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Overview  over  Discrete  Time  Domain  Methods  in 
Electromagnetic  Field  Computation 

Peter  Russer  * 

Abstract 

The  modelling  of  fields  in  the  time  domain  describes  the  evolution  of  physical 
quantities  in  a  natural  way.  Transient  phenomena,  nonbnear  and  dispersive  beha¬ 
viour,  the  characteristics  of  systems  with  moving  boundaries  or  with  time  dependent 
properties  are  best  described  in  the  time  domain.  In  this  contribution  different  ap¬ 
proaches  for  time  domain  modelling  of  electromagnetic  fields  are  compared. 

1  Introduction 

For  electromagnetic  field  modelling  numerous  techniques  have  been  developed  (1 ,2,3).  The 
modelling  of  fields  and  networks  in  the  time  domain  is  highly  attractive  since  it  describes 
the  evolution  of  physical  quantities  in  a  natural  way.  Time  domain  modelling  is  especially 
advantageous  m  the  case  of  transient  electromagnetic  fields,  fields  m  nonlinear,  dispersive 
or  time^dependend  media  or  in  regions  with  moving  boundaries.  One  of  the  main  advan¬ 
tages  of  time^domain  modelling  of  electromagnetic  fields  is  the  local  dependence  of  the 
field  variables  on  space  as  well  as  on  time.  Within  discretized  space  and  time  the  state 
of  the  field  in  a  given  point  and  at  a  given  time  depends  only  on  the  field  states  of  the 
neighbouring  points  at  previous  times  This  allows  a  highly  parallel  computation  of  the 
time  evolution  of  the  discretized  field. 

In  modelling  of  high  frequency  circuits  we  have  to  deal  with  the  network  as  well  with  the 
field  concept  (Table  1).  Whereas  the  field  has  a  spatial  structure  the  network  structure 
IS  topological.  However  if  the  field  is  described  by  a  discrete  set  of  base  functions  as  it  is 
done  for  example  in  the  method  of  moments  [4,5]  or  if  the  field  is  discretized  with  respect 
to  space  we  obtain  topological  relations  between  the  state  variables  of  the  field  This 
allows  to  apply  netwoik-lheor-tical  methods  to  field  problems. 


^Lehrstuhl  fur  Hochfrequenztechmk,  Technischc  Universitat  Munchen,  Arcisstrasse  21, 
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Table  1:  Concepts  of  Field  Theoiy  and  Network  Theory 


NETWORK 


FIELD 


Topological  structure 
Time,  Amplitude 


Spatial  structure 
Time,  Amplitude,  Space 


Continuous 


•  Analog  Network 


•  Electromagnetic  Field 


Discrete 


•  Digital  network 

•  Discrete  modelling 
of  analog  networks 


•  Cellular  Automata 

•  Discrete  Modelling 
of  fields 


In  the  following  we  shall  focus  our  attention  on  four  approaches  for  time  domain  modelling 
of  electromagnetic  fields. 


•  The  finite-difference  time-domain  (FDTD)  method 
e  The  transmission  line  matrix  (TIM)  method 

•  The  field  modelling  by  cellular  automata 

t  The  field  modelling  by  multi'dimensional  wave  digital  filters  (MDWDF) 


2  The  Finite-Difference  Time-Donaain  Method 

The  finite-difference  time  domain  (FDTD)  method  is  the  mathematical  approach  for  the 
solution  of  partial  differential  equations  (6).  The  partial  derivatives  are  simply  replaced 
by  finite  differences.  In  1966  Yee  has  first  given  a  finite-difference  time-domain  scheme 
for  solution  of  the  Maxwell  equations  (7,8,9).  In  the  FDTD  method  space  and  time  are 
discretized  with  increments  Af  and  At,  respectively.  The  field  component  placement  in 
the  FDTD  unit  ceil  is  shown  in  Fig.  1.  The  side  length  of  a  unit  cell  in  our  notation  is 
2AI 

Space  and  time  coordinates  are  given  by  x  =  /  Af,  y  =  mAI,  z  =  n  AI  and  t^k  At,  The 
FDTD  scheme  for  the  solution  of  the  Maxwell's  equations  is  then  given  by 
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Figure  1;  Field  components  in  the  FDTD  unit  cell. 
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+7  (( +  l,m,n  + 1)  -  (( -  l,m,n  +  1)  + 

(l,m  -  l,n  +  1)  -  (;,m  +  l,n  +  1)]  (6) 

with  the  stability  factor 

2cAt 
Af  ’■ 

where  c  is  the  velocity  of  light,  At  is  the  time  interv^,  A/  is  the  length  interval.  The 
condition  for  stability  is  given  by 

o<4  (8) 


The  FDTD  scheme  gives  the  new  field  state  at  the  time  at  as  a  function  of  the  field 
states  at  (k  ->  i)  At  and  (k  —  2)  At.  Also  the  new  spati^d  components  at  /,  m  and  n  are 
related  to  the  spatial  components  at/±l,/±2,  m±lym±2,  nd:l  and  n  d:  2.  However 
all  electrical  field  components  with  even  values  of  k,  are  only  related  to  the  magnetic  field 
components  with  odd  values  of  k  and  vice  versa. 

investigating  planar  circuits  within  the  magnetic  wall  mode)  a  two-dimensional  finite 
difference  scheme  may  be  applied  (10).  For  the  circuit  plam.  parallel  to  the  x  -  y-plane 
the  electric  field  exhibits  only  the  ^-component  and  the  magnetic  field  exhibits  only  the 
*-  and  y-  components  The  surface  current  J  flowing  in  the  upper  plane  of  the  planar 
circuit  is  given  by 

Ja-e,xH  (9) 

where  Ci  is  the  unit  vector  in  z  direction.  The  voltage  V  between  the  plates  is  given  by 


V  =  -<f£. 


where  d  is  the  distance  between  the  plates. 


(11) 

=  (12) 

Cl  IS  the  capacitance  and  Li  is  the  inductance  of  an  arbitrary  square  element  of  the 
two-dimensional  parallel  plate  line. 


+7-r  (V*(m  +  l,n)- V‘(»n  +  l,n)] 


j;^'(m,n)  =  J‘-'((m,n)  + 


+^(v‘(m,n  +  l)-V‘(m,n  +  l)] 


(14) 
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V‘+'(m,n)  =  V*-'((m,n)- 

-  ^  (4  (">  +  1. ")  -  4  <>”  -  1 ,  ti)  + 

+Jj(m,n  + 1)'- J‘(m,n  -  1)]  (15) 

Nonrectangular  grids  have  been  treated  for  the  PD  method  |U).  The  analysis  of  radiating 
structures  requires  the  termination  of  the  grid  with  ^sorbing  boundaries  [U,1243j' 

3  The  TLM  (Transmission  Line  Matrix)  Method 

The  transmission  line  matrix  (TLM)  method  was  developed  and  first  published  in  1971 
by  Johns  and  Beurle  (14,15,16,17).  In  the  TLM  method  the  physical  analogy  and  the 
isomorphism  in  the  mathematical  description  between  the  electromagnetic  field  and  a 
mesh  of  transmission  lines  is  used.  The  field  evolution  is  modelled  by  voltage  pulses 
propagating  on  the  mesh  lines  and  being  scattered  in  the  mesh  nodes.  Field  theoretically 
the  TLM  method  is  based  on  the  Huygens  principle  [IS, 19) 

In  the  TLM  model  space  and  time  are  discretized  in  length  intervals  Af  and  time  intervals 
At.  The  intervals  Af  and  Af  correspond  to  the  real  space  and  time  intervals  without 
scaling  if  Af  =  V2coAf/v/<7  is  chosen,  where  is  the  relative  permittivity.  Fig.  (2) 
shows  the  port  numbering  of  a  two-dimensional  TLM  shunt  node. 
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Figure  2:  A  two-dimensional  TLM  shunt  node. 

The  scattering  of  impulses  at  a  shunt  nodeis  described  by  the  following  equation: 


V, 

V, 

Vt 

=  s 

vs 

vs 

n 

v* 

k 

K 

(16) 
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Zo  is  the  characteristic  mesh  line  impedance,  given  by 

krjo 


Zo  = 


(17) 


where  rjo  =  3770,  h  is  the  height  of  the  planar  structure,  and  c,  its  permittivity. 
The  node  scattering  matrix  for  a  shunt  node  is 

I  1 


S  = 


-I  I 
a  a 
1  -1 
a  a 
1  I 
a  a 
I  1 


1 

a  a 
1  -1 


(18) 


Lossy  subregions  can  be  modelled  in  TLM  by  connecting  a  lumped  resistor  or  an  infinitely 
long  transmission  tine  stub  across  each  mesh  node  (20).  Also  the  modelling  of  nonlinear 
passive  elements  has  already  been  treated  [21,22,23].  Nonlinear  active  regions  may  be 
modelled  by  connecting  nonlinear  active  circuit  elements  to  the  mesh  nodes  (24,25) 

The  voltage  wave  puUes  incident  on  a  TLM  mesh  node  depend  on  the  voltage 

wave  pulses  fc-iV'Jw  „  emerging  from  the  neighboring  nodes  as  follows: 


Ki.m.n 

k+lKi, m.n  “k  ^.m.n 


on 


Eqs  (16)‘(18)  and  (19)  describe  the  complete  algorithm  for  the  time  discrete  Peld  evolu* 
tion. 

For  the  two-dimensional  case  Johns  has  shown  the  relation  between  the  FDTD  method 
and  the  TLM  method  [28] 

E,  =  pq’^  +  r 


with 


^inu) 


q 

-  (1*111 

r  =  -1 


*+i  =  9C  ^jfEt  +  r*V*^ 

£,  =  ijCp  „£,?CrCp  j.i  £,  +  «CrCp  .-i  V' 


(20) 

(21): 

(22) 

(23) 

(24) 

(25) 


With 


CrCr  =  al 


(26) 
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Figure  3:  Three-dimensional  symmetric  condensed  TLM  node, 
where  a  is  a  constant  we  obtain 

k+\Ei  =  qCp  kBtqCrCp  i.iE,ak-iEt  (27) 

W'.*  now  have  an  algorithm  relating  hi^’»  kE,  and  to^.i^,.  We  have  reduced  the 
'  ariables  but  increase  the  time  depths  of  the  algorithm 

Different  TLM  schemes  have  been  proposed  for  the  three*din’ensional  case  (16),  A  sym* 
metrical  three-dimensional  condensed  node  has  been  introduced  by  Johns  (26, 27],  Fig.  3 
shows  the  symmetric  condensed  TLM  node  In  the  three-dimensional  case  we  hava  to 
introduce  twelve  wave  amplitudes  The  voltage  wave  vector  is  given  by 

V'  =  [v;v;viv;viviv;v;viv{„v,\v,’,f  (28) 

the  incident  wave  pulses  *V'  at  t  =  and  the  scattered  wave  pulses  at  »  = 
(fc  +  l)Af  are  related  by 


H,V'  =  S.V 


(29) 


the  scattering  matrix  S  given  by 


0  T  T’' 
TT  0  T 
T  0 
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S  = 


with 


T  = 


0 

0  -i 
I  0 
i  0 


1  _1 


We  introduce  the  field  vector  F  given  by 


The  wave  amplitude  vector  fcVj  „  ^  is  related  to  the  field  vector  by 


wjth 


0 

0 

0 

0 

0 

0 

1 

1 

1 

1 

0 

0 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

0 

0 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

1 

0 

0 

0 

0 

1 

-1 

0 

0 

1 

-1 

0 

0 

0 

0 

-1 

1 

0 

0 

1 

1 

0 

0 

0 

0 

1 

-1 

0 

0 

0 

0 

(30) 


(31) 


(32) 


(33) 

(34) 


(35) 


Structures  to  be  analyzed  in  TLM  may  be  segmented  into  substructures.  This  method 
first  proposed  by  Kron  is  called  diakoptics  |29,30,31).  Within  diakoptics  the  scattering 
of  waves  by  boundaries  is  expressed  via  discrete  Green’s  functions  or  so-called  Johns 
matrices  [Hj.  Discrete  Green's  functions  may  also  be  computed  algebraically  (32)  The 
TLM  method  in  connedtion  with  absorbing  boundary  conditions  has  already  been  applied 
to  the  analysis  of  a  slot  antenna  [33] 


4  Cellular  Automata 


John  von  Neumann’s  most  extensive  work  in  the  theory  of  automa  .a  was  the  investigation 
of  cellular  automata  (34).  The  results  of  this  work  arc  contained  in  the  manuscript  “The 
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Theory  of  Automata*  Construction,  Reproduction,  Homogeneity”  (35).  John  von  Neu¬ 
mann’s  cellular  automata  are  homogeneous  two-dimensional  arrays  of  square  cells,  each 
containing  the  same  twentynine-state  finite  automaton.  Any  cell  can  assume  at  a  given 
time  the  unexcitable  state,  one  of  twenty  quiescent  states  or  one  of  eight  sensitized  states. 
The  unexcitable  state  represents  the  presence  of  no  neuron.  Quiescent  cells  can  respond 
to  stimuli  from  adjacent  cells.  The  state  of  a  cell  at  a  given  time  is  determined  by  a  set 
of  transition  rules.  The  state  after  the  transition  is  determined  by  the  initial  states  of 
the  cell  and  of  its  four  nondiagonal  neighbours.  John  von  Neumann  has  shown  that  the 
twentynme  states  are  sufficient  to  accomodate  all  logical  and  construction  circumstances 
that  may  arise  and  also  to  establish  all  transition  rules  for  moving  from  one  state  to  the 
other.  He  demonstrated  the  logical  universality  of  the  cellular  automata  by  showing  how 
l\iring’s  machine  model  can  be  reformulated  in  terms  of  cellular  automata. 

John  von  Neumann  also  has  planned  the  continuous  model  as  a  further  refinement.  In  1969 
Konrad  Zuse  discussed  the  modelling  of  the  Maxwell’s  equations  by  cellular  automata  [36]. 

Cellular  automata  now  meet  with  growing  interest  for  modelling  of  physical  systems  [38]. 
The  discrete  system  may  be  described  in  state  variable  form  (37).  An  autonomous  system 
is  described  by 

i{r,t+l)=  A{r',r)z{r',t)  (36) 


In  vector  notation  this  is  given  by 

2(f  +  l)  =  Az{0  (37) 

where  2(1)  is  the  state  vector  and  A  is  called  the  transition  matrix.  The  discrete  time 
variable  is  f,  and  r  is  the  o'serete  space  variable. 

As  an  example  we  consider  the  telegraphist’s  equation  [37j. 


The  corressponding  difference  equation  is  given  by 

-CA(V  s  <7v  +  A»i 


At  and  A,  are  given  by 


Aii(a:,f)  =  i(x,<  +  l)-i(x,() 
Ar*(x,f)  =  t(x  +  l.f)-t(x,() 
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The  state  vector  z  and  the  transition  matrix  A  m  eq.  (37)  are  given  by 


Using  the  Z  transformation  with  respect  to  x  the  variable  given  by 


(45) 

and  with 

vsO 

z(A.2(i)}  =  (r‘-i)zw*)} 

(46) 

we  obtain 

(47) 

with 

> 

II 

1 

1  •'«B 

(48) 

This  result  may  be  transformed  into  a  digital  circuit  Fig.  4.  shows  the  corresponding 
digital  model  of  the  lossy  transmission  line. 


5  Field  Modelling  by  Multi-Dimensional  Wave  Digital  Filters 


The  numerical  integration  of  partial  differential  equations  using  principles  of  multidi* 
mensional  wave  digital  filters  (MDWDFs)  has  been  proposed  by  Fettweis  [43,44|.  The 
continuous-domain  physical  system  is  simulated  by  means  of  a  discrete  paiaive  dynami¬ 
cal  system.  In  this  method  in  a  first  step  the  partial  differential  equation  is  modelled 
by  a  multidimensional  analog  circuit.  This  circuit  is  then  transformed  into  a  MDWDF 
equivalent  circuit  (45).  The  advantage  of  this  method  is  the  robustness  of  the  algorithm 
even  under  conditions  of  of  rounding  and  truncation  operations. 
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i(x,<)  t(i  +  l,<) 


Let  us  consider  as  an  example  again  the  two-dimensional  electromagnetic  field  problem 
associated  with  a  transverse  electric  field  between  two  metal  plates.  We  investigate  the 
equations 


dtt 


dh 

ih. 

di, 

di, 


(49) 

+  ^  =  “a(l) 

(50) 

(51) 

The  variable  h  corresponds  to  time,  whereas  <i  and  represent  the  spatial  coordinates. 
The  current  density  components  in  the  upper  metal  plate  in  the  directions  and  ti, 
respectively,  are  given  by  ij  and  tj,  and  v  is  the  voltage  between  the  metal  plates.  The 
inhomogeneous  terms  U|,  uj  and  U3  represent  distributed  impressed  sources. 

In  the  three-dimensional  complex  frequency  domain  we  obtain  the  algebraic  equations 


{psL  +  R)!i  •¥  PiHih  — 

(52) 

{p^L  4-  R)Ii  +  pzHih  = 

(53) 

Pi^sA  +  +  (paC  +  G)i^/3  =  U3 

(54) 

where  pi,  and  p3  are  the  complex  frequencies  The  corresponding  analog  circuit  repre¬ 
senting  the  partial  differential  equations  is  depicted  in  Fig,  5  Note  that  the  this  equivalent 
circuit  represents  the  complete  field. 


Figure  5'  Multidimensional  analog  circuit 
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using  the  wave  digital  filter  design  rules  the  circuit  according  to  Fig.  5  is  converted  into 
the  wave  digital  filter  circuit  shown  in  Fig.  6.  The  methods  of  MDWDFs  ensure  that 
the  algorithm  is  recursible  and  that  the  full  range  of  robustness  properties  of  WDFs  is 
conserved 


Figure  6:  Wave  digital  filter 


6  Conclusion 

We  have  compared  four  different  methods  of  discrete  time  domain  analysis  of  electroma* 
gnetic  fields.  The  methods  originate  from  different  theoretical  frameworks.  Whereas  the 
FDTD  method  is  based  on  mathematical  considerations,  TLM  originates  from  a  line  mo* 
del.  The  method  of  cellular  automata  is  based  on  the  theory  of  automata  and  systems  and 
the  method  of  MDWDFs  uses  the  analogy  to  multidimensional  circuits  and  wave  digital 
filters.  There  are  interesting  analogies  between  these  methods.  Each  of  these  methods 
contributes  special  advantages  and  interesting  contributions  to  the  solution  of  problems 
also  relevant  in  the  other  models 
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ABSTRACT 


"Ultrawldeband"  (UWB)  and  "Very  Short  Pulse"  (SP)  provide  alternative  designations 
for  the  same  class  of  transient  electromagnetic  wave  phenomena.  However,  UWB  relates 
these  phenomena  to  the  frequency  domain  whereas  SP  relates  the  same  phenomena  to  the 
time  domain  (TD).  By  generating  SP-TD  data  through  UWB  frequency  synthesis,  the 
evolution  of  the  TD  signal  with  increasing  bandwidth  can  be  traclcsd  systematically  to  its 
highly  localised  UWB  form.  Localized  pulse-lilce  feature:  (observables)  in  data  suggest  tliat 
modeling  and  i.nterpretation  in  terms  of  a  TD  "obsetvable>based  parametrizatlon”  (OBP)  is 
physieally  more  appropriate.  Implementing  a  TD-OBP  requires  new  thinking  and  concepts 
directly  in  the  time  domain.  A  systematic  OBP  strategy  for  learning  to  "think  TD"  via 
identification  of  TD  wave  objects  is  proposed  and  illustrated  by  examples  involving  SP 
excitation  of  layered  media,  strip  gratings  and  aperture^oupled  cavities.  Of  special  interest 
is  the  TD  evolution  of  the  strongly  dispersive  leaky  modes,  Floquet  modes  and  cavity  modes, 
and  their  role  in  synthesizing  the  SP  response. 
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source>exciied  pulsed  transient  fields  are  expressed 
direetfy  in  the  tune  domain  (w  spatial  spectral 
auper^Uon  of  rrmrienr  basts  neidi. 

2.  The  /oMd  wavefrcm-resonance  algorithm,  whldi 
constructs  t!me*d^>cndent  fields  scattered  by  targets 
or  inhomogeneous  media  in  terms  of  seV<wlstent 
comblnerions  of  proptSSlM 

oscOloMt  (fesonanee)  bssis  /»«!*•,  ^ 

svttematoed  the  fundamental  distinction  be^en 
the  earW  time  and  late  time  response,  and 
thereby  clarified  the  limitations  associated  with  the 
sin^irtcy  expansion  method 

3.  The  pKuse-spact  beam  algorithm  for  stlf<onslsteni 
^conposiuon  and  subsequent  recomblnauon  of 
time-hamoftic  or  transient  fields  into  wmdowed 
bean-type  basU  fields  on  I  eoafifuratioa* 
wavenumber  phase  space  lattice  or  ui  a  phase  space 
continuum.  The  oeam-type  fields  wa  good 
pr^gators*  »d  are  tracked  conveniently  through 
cempAx  propagation  and  satterlng  eneounters. 

4.  The  complex  sowee  pains  (CSP)  method  for 
modeling  of  nuired  bean  type  inputs,  and  the 
scattenhg  ei  these  pulsed  beams  by  eovtronmems. 
The  OBP  modeling  strategy  described  above  has  been 

applied  to  a  vancty  of  propagatioo  and  scattenng 

esvironments  in  electromagnetics  and  underwater 

_ .!„  f  'rvacA  aMMtiMflnn*  Vawfi  h«*n  nniT’&nlv 


subsequent  arect  iti'vro,  »  piw«»«7  wv*-# 
i^emcnied  for  various  environmental  cowuons.  Sr 
plane  wave  seaitertng  from  a  sl‘t  coupled  cavity  has 
already  been  implemented  (W].  Other  examples  include 
SP  adutioa  of 

a.  tingle  and  multi^e  strip  targets 

b.  periodic  and  quasi  periodic  arrays  of  stiipe  without 
and  with  a  rcQMting  plane  backup. 

e.  planar  and  cyllndrically  layered  dielectrics 

These  Utter  examples  pennit  investigation  of  the  TD 
buildw  of  multiple  inteiacuons  and  of  disptf  uve  effeets 
associated  with  periodic  structures  and  guided  motto: 
underttafiding  tms  phenomenoloty  b  expected  to  be 
useful  for  subsequent  studies  of  Sr  scattenng  by  clutter 
and  W  to  clutter.  Preliminajy  resulu  are 

presented  here. 

2.  SytithuU  opilena  for  TD  dyadic  Gma’i  Funetlons 
A.  Plana  layered  dielectrle  madia 
i.  Fomal  aspects 

Dyadic  Oreen'i  fooctions  for  plane  stratified 
diek^  may  tend  to  perfto  coa^toia  ia  tte 


funedoBS,  whidi  thereby  play  a  fundamental  role  m  field 
synthesis.  With  *  and  ;  ■  (it,y)  denoting  the  coordinates 
perpoidicular  and  parallel  to  the  hierfoces, 
respectively,  the  potemisJ  sohitiMS  may  be  constructed 
by  spectral  (Fou^)  decompositioA  into  the  transverse 
wSmnumber  domain  l^t  •  (li.lv)  ^  temporal 
wavenumber  (frequency)  domain  ca  s^g  the  z* 
de^dent  problem  in  the  ( k ,  .w)  domain  by  invoking  the 
boundary  conditions  aertts  layer  interfaces  and  the  imtiil 
conditions  at  the  source  loatlon  r '  ■  (x'.y',z*),  and 
thereat  performing  spewil  synthesis.  This  to 
the  foilowing  gen^c  form  for  the  potential  function  f  at 
the  space*time  obsenraiion  point  (;  ,t), 

'.«')•  ^Fll  <ii  •'(f  .1*  M'  w !? ,)  (1) 

where  the  wave  o^cct  F,  denited  in  the  (spacc*iime)' 
(wavenumber-frequenQf)  phase  space,  has  tbe 
condition 

(la) 

and  tbe  denominator  D  is  written  conveniently  as 

(H>) 

Here,  R,  and  R.  are  spectral  domain  reflation 
cocliclents  Seen  looldng  along  the  (+z)  and  (‘Z) 
directiotu,  respectively,  from  the  source  plant  z  •  z  . 
For  the  most  venatHe  treatment  of  the  spectral  integrals 
m  0).  it  is  useful  to  separate  (1/D)  lelf.conststently  into 
hytad  *fay‘mode*  constttuenu 

iy>.  S  (R.R.)“.l)-(R.R.)’''+(R.!<.)’'’‘’P‘. 

a*H> 

|R.Ril<|.  (2) 

Obtained  by  partial  power  series  expansion,  each  term  in 
the  series  portion  resprcsenis  n-times  mfiected  waves 
whleh  eanV  maniputaied  into  generalized  ray  fields 
after  synthesis,  whereas  the  remainder,  with  iu  resonant 
denondnalor,  can  be  manipulated  into  guided  mode 
fields.  Note  that  the  amplitude  of  the  modal 
depends  self*coniisieBtly  on  tbe  indexes  N;  and  Nt  of  the 
l^ed  rays,  with  Nj  ■  0,  Nj  ■  »  hrmshlng  a  pure  ray 
field  series,  and  iniact  retention  of  O'*  (urrushing  a  pure 
mode  field  series. 


formal  i, . , 

ueat  the  real  potential  f  — 
in  tbe  kwer  half  of  the  complex  i  plane. 

•  • 

TUm)  -  7  f  ♦  (..»)  e^  dt;  f(.  t)  -  Re  f*  (.1)  (3a) 

•o» 

In  the  reduction,  the  most  significant  options  are 
assod^  with  v^tber  the  aeinversion  is  performed 
i^iir  or  he/bre  the  b  inversioa  The  former  is  the 
tomentitmd  route,  going  from^the  fuU  solution  in  the 
frequency  domain  to  the  tune  domain,  while  the  latter 
i:  jfo  is  noncomendvtal  and  utilizes  transient  plane  vaves 
la  the  point  source  lyniheaU.  Exercisbg  these 
alternaiivM  together  with  the  decomposition  in  (2),  and 
p^onning  contour  deformations 
pknes  ioste^t  descent  paths  (SDP)  through »Mle 
DOW  of  the  mtegrands.  one  may  arrive  at  exaa  hybrid 
r^>roode  formulaiions  eompnseo  of 

Urn  ^  (generalized  r^nelds)« 

4  *2  '(modified mode  fieMi)^ 
qCNu) 


J 
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*  ({Uy>mode  remtinder)N^  (4) 

Etch  n*term  m  the  first  series  contiins  in  SOP  iQtetnJ 
Mth  a  *ray  phase':  its  a^ptotie  evaluation  by  the  saddle 
'4nt  method  yields  a  space*tnne  field  (wavefront), 
'^ecb  o-term  in  the  second  series  arises  from  a  residue  at 
'le  pcto  of  D*'  encountered  donnc  the  SOP  comour 
oeforxnatioiv  whQe  the  last  term  in  (4)  is  a  ^brid  spectral 
tmecral  remainder  alone  the  N]  and  N]  SOP'S,  ha^e  a 
ray  saddle  point)  phase  but  a  moik-t^  amplituM. 

Both  the  second  and  third  terms  in  (4)  are  derived  from 
the  second  term  in  (2).  Details  of  lo^  reductions  may 
be  found  in  (32).  ' 

Ibe  conventional  and  nonconvendonal  routes  In 
performini  the  spectral  ^thesis  lead  to  the  foDowfrie 
alternative  treatments  of  the  dispenion  etjuation  tbit 
defines  the  modal  contn'butions  In  (4): 
•)D<Hl!,)*0->!F|^Hk)D(u„!!_).0»<.w(lC|)  (5) 

^e,  ^  q*p  idratifiea  the  conventioita]  fre^ency 
'v),  whereas 


IS 


This  relation  can  be  schematized  by  the  $paee>u'me 
dependent  self*conai«ent  ray»mode  trajectories  to  Re.  2, 
Nwih  V,  rmesentini  the  modal  pmp  speed.  The  rnodal 
phase  evamaied  at  ^i.ku)  becomes 

where  v^  the  modalphore  speed 

X  Mtfnencal/eru/(r 

Numerical  resuUi  for  the  eoofifuratien  In  Fit.  2  have 


depewcst  spatial  wavenumber  p^es  k  (w).' 

>m  idcmines  the  nonconvention^*^  wavenumber 
tpendent  frequency  poles  wm  0;,}.  The  eorrespondlni 
m^  phases  and  amplitudes  may  be  inferred  from  the 
^ral  teietrands  k  (It),  with  (2).  evatuaied  at 
(w.k^(«))ana(i.^(kj),k^).  respecuveV 

i  Modt  asymptones 

The  spectral  intcirsls,  which  remain  in  the  complex  u 
and  the  complex  h,  plencs  for  the  options  in  <5a}  and 
(SbX  respectively,  can  be  evaluated  asymptotically  ^  the 
saddle  point  method  The  rcspccthe  sa^U  point 
conditions  w«(i^andl^^«le^^are  specified  by 

a)<»,(.u)/4.|*  .0;  ♦,(•!:,) I,.«0  («) 

»b«e  <«iioia  Uie  pidicai  opetMOt  In  tlK  ); 
Spectral  domaia  This  implies  that  conditions  (5)  and  (6) 
must  be  aatisfled  slmuhmtousfy  to  yteld  the  fiwd  space* 
tune  modal  spectral  wavenumbers  w«(j.t7',t*)  and 
lU  'll');  ihese/mwf  values  are  Independent  of  the 
maoner  (cooveotlonal  or  nonconvemional)  by  which  they 
were  derived.  The  simultaneous  requirements  ib^  can 
be  sebematized  by  fraphical  construction  whj^  locates 
on  the  four*diraensiooal  dispersion  surface  0(i.i.k)»0 
those  points  (w^k  ^)  that  have  a  surface  normal  p  pv allel 
to  the  fouT'dincnsiona]  space’iime  riy  vector 
(T'I  .H  X the  oriemationi  of  the  eoer^te  axes 
la  the  spectral  and  confi|urai!onal  domains  arraneed  as 
m  Fig,  \  (33},  Here,  k  “(ijikjCk^))  is  the  three- 
dimcoiioDal  waYcvcctur.  WhSc  )  are  generally 
cam;^  go^  ^oximitions  for  the  phases  of  weakly 
complex  tnodes  can  be  obtained  ^  plotting  the  real  parts 
(c^.k^)aitnrjg.J, 

^  Consideratioa  from  here  on  wiU  be  restricted  to  a 
s^e  diejecuK  layer  on  a  perfectly  conducting  ground 
plaac.  wstb  the  source  pomt  (r  *,t*)  ■  (0,0^a‘,0)  umd^ 
and  the  diceivation  point  (r,t)  located  outside,  the  layer 
(R&  2X  The  relevant  mode  fiel^  arc  now  the  leab 
audes  beauM  ite  trapped  »od*s  have  eiiernd 
fields.  ^  saMe  pouil  conditions  w  (d), 
ehich  yield  Idenucal  residts  f(K  mo^  q«p  or  m,  are 
ei^UdtJy 

tttU  U  U 


(7) 


d^eet  Summation  of  n^fele  reflected  wavefronts  (Hrst 
term  in  (4)  with  Ni  •  a  Ki  -  N)  and  b)  via  leeb  inode 
summatiM  (second  term  m  (4)  with  the  full  spectral 
in  (la))  over  rnooet  whosa  frequencies  ^ 

|e  Mtbin  tbe  frtqueaQr  window  of  the  pulse  speeirum 
F(w),  BeaiM  the  multiple  reflected  wivefroi)iTiel&  in 
a)  ara  nondlspertiva,  they  can  be  evaluaied  4K  simple 
fecTo;  tbe  wavefront  scries  is  truncated  direetty  W  the  n* 
flayed  arriviil  times  that  Ht  Into  the  tune  tnietval  ^ 
obMrvatioa  andjnduectly  by  the  d^ease,  at  each  air* 
d»lecinc  reflection,  of  the  field  amplitude.  The  results 
ara  shown  in  Figs.  ^  aM  are  npluned  in  the  figure 
espitets.  Of  special  interest  it  the  synthesis  of  the 
.  wavelront  dnvan  $P  observables  by  successive  addition  ^ 
tix  tmiai/  dttpwuvt  lulw  motel  tliij  Uluiinm  ibt 
role  of  dupersion  under  Sf  coiwiitioni  iiul  the  vutue  of 
OBP  in  the  fuLV  lyiubesiz^  data. 

The  above  analytic  summary  and  numerical  samples 
have  been  extracted  from  a  full  manusaipt  bcmi 
preparedfor  publieation  [32). 
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Fig.  2.  Leaky  node  dispenioa  sut^m  for  c^aoe  leered 
didectric,  tad  craphlcid  construeih'D  for  determination 
of  the  rM  modal  wavenumber  k,t.i  .t)  and  fre^iefiQr 
<^((f  lO  defined  in  (6).  Ihe  subrenpi  t  been  omitted 
and^iS}:,. 


1  Korizofttal•eleetri^dipole>ex€tted  dieleeuic  laver 
with  relative  peronttM^  ci  relative  pemeawliQr 

bai^  to  a  peneoV'^ndoctiPf  (FEC)  ground 
pt^  with  owetver  located  in  the  outside  vacuum. 
Schema tinttion  of  aelf>co«tstem  ray^mode  traiectones 
that  eitiUiah  the  time^donain  asyxnpto^  kato  mode 
Held,  as  nested  in  (^.  Both  the  Inodent  neid  and 
the  detacniag  leaky  ray  field  to  the  observer  are  phase 
matched  to  the  leaW  mode  loogttodinal  wavenumber; 
tfau  matehiAf  eondltioa  deiemmes  the  ^ee-time 
dependent  reactive  ray  iDglei  9ii(r  ,t)  and  9i(r  >t)* 


Fig.  1  Numerical  results  for  TE  coninbut/on  to  ^  field  component,  normalised  to 
incident  pulse  ampliiudc  E«.  Physical  conCguraiioo  as  in  F:g  2.  with  problem 
parameters  listed  on  each  figure. 
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a)  i(^  pulse  [analyiic  Reyieigh  wavelet  ft(r)«it| 
and  iu  fremieocy  spectrum  (F(v)«(a<r 
ti«(XSB/cj.  flt)«Ref*(rX  the  circles  on 
the  uB  tsaB/c  aids  tdrati^  the  cutoff  Ire^cnaei  of 

the  l^mooes,.Oaty  the  first  22  modes  (0^m$ll)are 
covered  by  the  pulse  spectrum,  with  m»3  excited  most 
strong^. 


Direct  vavefrent  caleuUtieA 


Ex/EO  Ex/EO 


LB.  Febftn.  V«y  Short  Puhe  Satteriat:  Tune  Domain  Ptrametrization 


7  /r«m  msl  to  ms?  &  modre.  from  mel  to  msS 


«)  buildgp  of  received  ii^nal  by  mod&i  ti^rposilioiv  mode  m«3  (most  ttron^Ijr 
excited):  modet  mc«e$  l^m<7:  model  with  detult  ck 

eoastructtve  end  destructive  iaterferenee.  Mode  sum  O^mSU  compleieV 
recoAstrucu  the  Held  in  b).  Note  the  deiced  tum>on  tunes  »t  the  observer;  thus  are 
compatible  with  the  earliest  phase  matched  propafatim)  paths  identified  io  Fi^.  2. 
The  abrupt  start  should  not  be  taken  Kterally;  It  ames  from  nonmiform  asymptoto. 


Finite-Difference  Time  Domain  Moclellino 
of  Eleotromaonetlo  Fleicis 

Inoo  Wolff.  Oulfbura  Univtnitv 


for  (lie  dosign  of  planar  microwave  integrated  circuits  up  to  1985  mainly  analysis 
techniques  m  the  frequency  domain  have  been  used  With  the  requirements  for 
new  and  flexible  tools  in  the  design  of  planar  circuits  eg  with  closely  coupled 
element*)  and  three'dimensional  diicontmuitioc  iiko  oirbridgco.  alternative  techniques 
must  be  studied.  One  of  these  techniqt.6s  1$  the  finite-difference  time  domem 
analysts  (FOTD)  which  In  principle  1$  known  since  30  years.  Yee  already  in  1966 
proposed  this  technique  for  the  analysis  of  electromagnetic  boundary  value  problems 
(13 

During  a  long  time  the  fOTD  technique  only  was  used  to  quantitatively  demonstrate 
electromagnmin  field  .<inlitt(on  in  th^  time  domain  Only  tha  introduotion  of  oooorbmg 
walls  made  this  technique  to  a  powerful  tool  for  the  solution  of  real  problems. 

The  FDTO  is  a  numerical  method  for  the  solution  of  electromagnetic  field  problems 
which  has  a  large  numerical,  but  a  low  analytioal  expensa,  Daapite  the  large 
numerical  expense  it  is  believed  to  be  one  of  the  most  efficient  techniques,  because 
basically  it  stores  only  the  field  distribution  at  one  moment  in  memory  instead  of 
working  with  a  large  equation  system  matrix.  The  field  solution  for  each  other 
time  then  is  determined  from  Maxwell's  equations  and  is  caloulated  using  a  time 
stepping  procedure  based  on  the  finite-difference  method  in  time  domain.  Available 
algorithm,  called  the  'leapfrog  algorithm*  fits  very  well  on  modern  computer  archi¬ 
tectures.  so  that  the  data  required  to  describe  e  three-dimensional  field  distribution 
can  be  handled  m  a  reasonable  time.  Therefore  It  can  efficiently  be  implemented 
on  vector  or  parallel  computers  as  well  Sufficiently  accurate  results  can  be  received 
by  using  a  single  precision  floating  point  expression  requiring  only  four  bytes.  It  is 
a  further  advantage  that  the  transient  analysis  delivers  a  broad  band  frequency 
response  in  one  single  computation  run. 

In  this  talk  it  shall  be  demonstrated  that  the  technique  can  be  applied  to 
real  microwave  circuit  design  problems  It  will  be  shown  that  this  technique  enables 
to  mode)  arbitrarily  shaped  planar  structures  with  multiple  coupled  discontinuities 
and  planar  lines  and  three-dimensional  circuit  structures  Several  applications  to 
realistic  problems  of  modern  monolithic  Integrated  microwave  circuit  design  problems 
will  be  demonstrated  In  the  future  the  application  of  FDTD  method  surely  Is  a 
powerful  analysis  technique  for  nonlinear  microwave  integrated  circuit  design  by 
combining  physical  semiconductor  models  whtch  work  in  the  time  domain  with  the 
FDTD  description  of  the  passive  circuit  elements. 

(ll  KS  Yee,  '’Numerical  solution  of  imtial  boundery  value  problems  involving 
Maxwell's  equations  m  isotropic  metfta*.  JBEB  Trans  Antennas  Fropagat.  Vol 
U,  1966,  pp.  302-307 
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ABSTRACT 

In  this  workshop  paper  the  principal  features  of  TLM  analysis  of  electromagnetic 
fields  will  be  summarized,  and  research  trends  in  this  area  will  be  discussed.  Time  domain 
modeling  in  general,  and  TLM  modeling  in  particular,  is  focusing  on  the  realization  of 
a  new  generation  of  time  domun  simulation  tools  which  link  geometry,  layout,  physical 
and  processing  parameters  of  a  microwave  or  high  speed  digital  circuit  with  its  system 
specifications  and  the  desired  time  and  frequency  performance,  including  electromagnetic 
susceptibility  and  emissions.  Such  CAD  systems  will  most  likely  employ  dedicated  parallel 
processors  configured  in  a  3D  array.  Rirthermore,  the  specific  nature  of  discrete  time 
domain  algorithms  affords  optimization  and  synthesis  procedures  which  differ  radically 
from  those  employed  in  traditional  frequency  domain  CAD  tools. 

X  PROPERTIES  OF  TIME  DOMAIN  FIELD  MODELS 
1.1  Time-Stepping  Algorithms 

Most  time  domain  field  models  describe  only  the  local  properties  of  the  propagation 
space.  The  most  current  forms  are  based  either  on  a  discretization  of  Maxwell’s  Equations 
(Finite  Difference  -  Time  Domain  or  FD«TD  formulation)  or  on  the  description  of  space  by 
a  discrete  spatial  network  ('IVansxnission  Line  Matrix  or  TLM  formulation).  Fig.  I  shows 
the  basic  2D  TLM  impulse  scattering  process  which  can  be  considered  as  a  computer 
implementation  of  Huygens’  principle.  Finite  Element  formulations  m  the  time  domain 
are  alv>  possible  but  have  not  been  used  extensively  so  far. 

FD*TD  and  TLM  methods  employ  similar  but  different  formulations.  While  FD-TD 
is  expressed  in  terms  of  total  electric  and  magnetic  field  conponents,  TLM  uses  incident 
and  refiectfd  wave  quantities  in  a  spatial  network.  As  a  general  rule,  both  formulations 
are  equivalent;  for  each  TLM  scheme  there  eidsts  an  equivalent  FD'TD  formulation.  Fig. 
2  shows  two  such  pairs.  Figs.  2a  and  b  show  Johns'  distributed  TLM  node  |1]  and  Yee’s 
unit  FD*TD  cell.  |2].  Figs.  2c  and  d  compare  Johns’  condensed  TLM  node  [3)  and 
the  equivalent  FD*TD  scheme  derived  by  Chen  et  al.  (4).  Fig.  3  shows  the  dispersion 


diai&cteristics  the  discretization  schemes  in  Fig.  2  as  derived  by  Nielsen  and  Hoefer  [5J. 
For  low  frequencies  the  dispersion  surfaces  form  a  unit  sphere  in  all  cases.  However,  at 
higher  frequendes  the  dispersion  characteristics  the  condensed  TLM  node  and  Chen’s 
FD-TD  scheme  (Fig.  3b)  are  far  superior  to  the  other  two  (Pig.  3a). 

Clearly,  the  equivalent  TLM  and  FD-TD  schemes  possess  identical  dispersion  and  er¬ 
ror  characteristics  They  can  also  be  derived  formally  one  from  the  other.  Rirthennore, 
optimized  codes  for  eqmvalent  schemes  require  a  similar  computer  memory  and  execution 
time.  Nevertheless,  they  have  their  respective  advantages  and  disadvantages  when  im* 
plementing  boundaries,  disperdve  constitutive  parameters,  and  nonlinear  devices.  In  the 
final  analysis,  the  choice  between  TLM  or  FD-TD  is  based  mere  on  personal  preferences 
and  familiarity  with  one  or  the  other  method  rather  than  on  objective  criteria.  In  the 
following,  the  salient  features  of  time  domain  simulators  will  thus  be  described  in  terms  of 
TLM  formalism  with  the  understanding  that  there  exists,  or  could  be  found,  an  equivalent 
FD'TD  formulation  unless  indicated  otherwise. 

1.2  Requirements  for  Time  Domain  Field  A  nalysis 

The  principal  advantages  of  modeling  electromagnetic  fields  in  the  time  domain  are 
well  known.  However,  in  order  to  exploit  them  in  a  practical  application,  dispersive  and 
nonlinear  properties,  moving  boundaries,  and  sophisticated  signal  processing  procedures 
must  be  implemented,  which  include  forward  and  inverse  Fourier  transform,  convolution, 
and  absorbing  boundaries.  Another  important  requirement  for  practical  applicabihty  is  a 
graphic  user  interface  for  3D  geometry  editing,  parameter  extraction  and  display,  as  well 
as  dynanuc  visualization  of  fields,  charges  and  currents. 

The  feasibility  of  these  features  has  been  demonstrated  both  m  TLM  and  FD*TD  en> 
vironments  [6]-[8].  However,  the  computational  nquirements  for  irodeling  complex  struc- 
tures  with  such  methods  are  still  extremely  severe.  Therefore,  research  efforts  are  being 
focused  on  the  development  of  accelerating  techniques,  the  most  important  of  which  will 
be  <h$cussed  below, 

2  ACCELERATING  TECHNIQUES  IN  TLM  MODELLING 

In  the  following,  the  most  important  accelerating  techniques  will  be  briefly  described. 
The  first  exploits  the  localised  nature  of  the  time  domain  algorithms  through  parallel 
processing,  the  second  is  based  on  the  numerical  processing  of  the  time  domain  output 
signal  using  the  Prony-Pisarenko  Method,  and  the  third  involves  the  reduction  of  the 
so-called  coarseness  error  by  improving  the  properties  of  the  discrete  TLM  mesh  in  the 
vicinity  of  sharp  comers  and  edges. 


2.1  Parallel  Processing 

The  prindple  of  causality  ensures  that  any  diange  in  the  state  of  a  TLM  node  affects 
only  its  immediate  neighbours  at  the  next  computational  step.  This  allows  the  implemen¬ 
tation  of  TLM  in  a  form  quite  different  from  the  pr^am  on  a  serial  machine.  Since  in 
a  parallel  computer  each  processor  has  its  own  memory,  it  is  practical  to  assign  to  each 
of  them  an  impulse  scattering  matrix  and  a  set  of  boundary  c(mditions.  The  impulse 
scattering  matrix  incorporates  the  local  properties  of  the  computational  space  such  as 
permittivity,  permeability,  conductivity,  and  mesh  size  in  the  three  co-ordinate  directions. 
The  boundary  conditions  specify  whether  there  are  boxmdaries  between  a  node  and  its 
neighbours,  or  whether  the  nodes  are  connected  together.  This  parallel  implementation 
greatly  facilitates  variable  mesh  grading,  conformal  boundary  modeling,  and  the  simulation 
of  lughly  inhomogeneous  materials  and  complicated  geometries. 

Fig.  4  compares,  on  a  logarithmic  scale,  the  improvements  made  over  the  last  year  in 
computing  speed  using  various  programming  techniques  (9)  and  parallelisation.  The  orig¬ 
inal  matriz  formulation  by  Johns  [3]  requires  144  multiplications  and  126  additions  and 
subtractions  per  scattering  per  node.  Through  manipulation  of  the  highly  symmetrical 
impulse  ccattering  matrix,  Tong  and  Fujino  [9]  have  reduced  the  scattering  to  six  multi¬ 
plications,  66  additions/subtractions  and  12  divisions  by  four,  increasing  computing  speed 
over  six  times.  Programming  in  Assembler  rather  than  €4*4*  accelerates  the  process  again 
four  times,  Finally,  parallel  processing  increases  speed  by  more  than  two  orders  of  mag¬ 
nitude  over  the  fastest  serial  version  implemented  on  a  386  computer  in  04*4*  language. 
The  combined  measures  effectively  reduce  computation  times  from  hours  to  seconds. 

TIls  comparison  strongly  suggests  that  future  implementations  of  time  domain  simu¬ 
lators  for  CaD  purposes  will  be  based  on  dedicated  parallel  processors  or  supercomputers 
that  emulate  pnrdilel  processing. 

2.2  Signal  Processing 

The  fast  Fourier  LVansform  (FFT)  is  the  most  frequently  used  signal  processing 
methc**!  .or  extracting  the  spectral  characteristics  of  a  structure  from  a  time  domain  simu- 
Idtios. .  For  ciHcient  computations  it  is  of  prime  importance  to  reduce  the  number  of  time 
samples  req oiTed  to  extract  a  meamngful  frequency  response.  To  achieve  this  goal,  process¬ 
ing  of  the  lime  response  using  the  Pronj’-Pisarenko  method  has  been  applied  successfully. 
(10). 

In  this  appr.>ach  the  discrete  time  domain  output  signal  is  treated  as  a  deterministic 
signal  drowned  in  noise.  (Fig.  5).  The  signal  is  then  approximated  by  a  superposition 
of  dumped  exponential  functions  (Prony’s  method),  and  the  noise  is  minimized  using  Pis¬ 
arenko’s  model.  This  signal  processing  technique  reduces  the  number  of  required  time 
samples  by  typically  one  order  of  magnitude. 
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2.3  Reduction  of  Coarseness  Error 

One  of  the  principal  sources  of  error  in  the  TLM  analysis  of  structures  with  sharp 
edges  and  comers  is  the  so-called  coarseness  error.  It  is  due  to  the  insufficient  resolution 
of  the  edge  fidd  hy  the  discrete  TLM  network.  The  error  is  particularly  severe  when 
boundaries  and  thdr  comers  are  placed  halfway  between  nodes  as  shown  in  Fig.  6.  It 
is  clearly  seen  that  the  nodes  situated  diagonally  in  hront  of  an  edge  are  not  interacting 
directly  with  the  boundary  but  receive  information  about  its  presence  «dy  across  their 
neighbours  who  have  one  branch  connected  to  it.  The  network  is  thus  not  sufficiently 
"stifT  at  the  edge,  and  results  obtained  are  always  shifted  towards  lower  frequencies.  The 
classical  remedy  for  this  problem  is  to  use  a  finer  mesh  in  the  vicinity  of  the  edge,  but 
this  introduces  additional  complications  and  computational  requirements.  On  the  other 
hand,  the  dispersion  characteristics  of  the  condensed  3D  TLM  node  (see  Fig.  3b)  are  so 
good  that  the  velocity  error  is  practically  negh^ble  even  for  rather  coarse  meshes.  A  much 
better  and  more  efficient  way  is  thus  to  modify  the  comer  node  such  that  it  can  interact 
directly  with  the  comer  through  an  additt<mal  stub  as  shown  in  Fig.  7  for  the  2D  case. 
Smce  tlus  stub  is  longer  than  the  other  brandies  by  a  factor  >/2  it  is  simply  assumed  to 
have  a  correspondingly  larger  propagation  velocity.  In  the  3D  case  up  to  three  stubs  must 
be  added  depending  on  the  nature  of  the  comer.  The  effect  of  this  comer  correction  is 
demonstrated  in  Fig.  8  which  shows  typical  results  for  the  first  resonant  frequency  of  a 
cavity  containing  a  sharp  edge  as  a  function  of  the  mesh  parameter  The  parameter 
p  is  proportional  to  the  fraction  of  power  carried  by  the  fifth  branch  of  the  comer  node 
and  is  equal  to  half  the  characteristic  admittance  of  the  comer  branch  when  normalized 
to  the  link  line  admittance  (see  Fig  7).  For  p  »  0  (no  comer  correction)  the  coarseness 
error  increases  almost  linearly  with  increasing  A/,  while  for  p  =  0.1  the  frequency  remains 
accurate  even  for  a  very  coarse  mesh. 

3  BOUNDARIES  IN  ARBITRARY  POSITIONS 

3.x  Accurate  Dimensioning  and  Curved  Boundaries 

The  accurate  modeling  of  waveguide  components,  discontinuities  and  junctions  re¬ 
quires  a  precision  in  the  positiomng  of  boimdaries  that  is  identical  to,  or  better  than  the 
manufacturing  tolerances.  If  boundaries  can  only  be  introduced  either  across  nodes  or 
halfway  between  nodes,  then  the  mesh  parameter  A/  would  have  to  be  very  small  indeed, 
leading  to  unacceptable  computational  requirements.  Similar  considerations  apply  when 
curved  boimdaries  with  very  small  radii  of  curvature  must  be  modeled.  It  is  therefore 
important  to  prori!de  for  arbitrary  positioning  of  boundaries.  The  basis  for  this  feature 
has  been  described  already  in  1973  by  Johns  (11). 
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Fig.  9  shows  the  concept  arbitreiy  wall  positioning  in  two'dimensional  TLM.  The 
boundary  branch  which  has  a  length  different  {torn  bXfl  is  simply  replaced  by  an  equivalent 
branch  of  length  A//2  having  the  same  input  admittance.  This  ensures  synchrcmism, 
but  requires  a  different  characteristic  admittance  for  the  boundary  branch  and  hence,  a 
modification  the  impulse  scattering  matrix  of  the  boundary  node,  (see  [U]).  The 
effect  of  such  boundary  tuning  is  shown  in  Fig.  10  which  indicates  that  the  length  of  the 
boundary  branch  can  be  continuously  tuned  over  a  range  of  more  than  one  mesh  parameter 
length  At  without  appreciable  error.  This  important  technique  removes  the  restriction  that 
dimensions  of  TLM  models  can  only  be  integer  multiples  of  the  mesh  parameter. 

An  alternative  method  which  avoids  the  modification  of  the  S-Matrix  of  the  boundary 
nodes  is  to  replace  the  extension  of  a  boundary  beyond  its  standard  position  by  an  equiv¬ 
alent  reactance.  The  differential  equation  cf  that  reactance  is  discretized,  resulting  in  a 
recursive  formula  for  the  impulse  reflected  by  the  boundary.  This  method  is  preferrable  for 
a  serial  type  computer  implementation  while  the  former  is  more  appropriate  for  a  parallel 
version. 


3.2  Moving  Boundarie.^  ard  Time  Domain  Optimization 

Since  it  usually  takes  considerable  time  to  bmld  up  a  quasi-stationnary  field  in  a 
structure  of  high  Q-factor,  optimization  based  on  a  new  complete  analysis  after  every 
modification  is  extremely  time  consuming.  Instead,  techniques  for  continuously  vaiymg 
the  boundary  position  and  other  characteristics  of  a  structure  during  a  TLM  simulation 
will  be  developed.  Two  different  methods  will  be  investigated.  One  is  to  modify  the 
scattering  matrix  of  nodes  situated  close  to  a  boundary,  the  other  is  to  generate  the 
impulses  reflected  by  moving  boimdaries  uring  recursive  algorithms.  In  order  to  implement 
automatic  optimal  tuning  these  measures  will  be  coupled  with  appropriate  optimization 
strategies.  Furthermore,  if  optimization  criteria  arc  to  be  formulated  in  the  frequency 
domiun,  a  sliding  Fourier  transform  window  will  be  introduced  as  well  in  order  to  extract 
the  time- varying  frequency  domain  characteristics  from  the  evolving  time  domain  response. 

4  NUMERICAL  SYNTHESIS  BY  REVERSE  TLM 

It  has  been  shown  recently  by  Sorrentino  et  al.  fl2]  that  the  TLM  process  can  be 
reversed  without  modification  of  the  algorithm,  yielding  the  source  distribution  from  the 
resulting  field  by  going  backwards  in  time.  Direct  numerical  electromagnetic  synthesis  is 
completely  unchartet  territory  as  yet,  and  the  exact  procedure  and  its  implementation  are 
not  very  dear.  The  desired  charactenstic  of  a  structure  or  component  is  usually  given  for 
a  limited  frequency  range  and  for  the  dominant  mode  of  propagation.  This  information 
is  insufiicient  to  synthesise  the  exact  topology  of  the  structure.  Therefore,  the  missing 


information  must  be  generated  and  added  the  designer.  Implementation  will  most 
likely  be  an  alternate  sequence  of  analyses  and  syntheses  which  will  converge  much  faster 
than  repeated  analysis  and  optimisation  in  traditional  CAD. 

5  CONCtUSION 

Computer  time  and  memory  required  to  model  realistic  electromagnetic  structures 
are  still  obstacles  when  it  comes  to  practical  applications  cf  time  domain  modelling  tech* 
niques.  Therefore,  considerable  research  eff^s  are  concentrating  on  ways  to  reduce  the 
computation  count  signihcantly.  In  this  workshop  paper  we  describe  three  di^erent  ways  to 
achieve  this,  namely  parallel  processing,  Prony*Pisarenko  signal  processing,  and  coarseness 
error  compensation  at  sharp  comers  and  edges.  All  these  methods  can  be  combined  to 
accelerate  TLM  simulations  by  several  orders  of  magnitude.  Since  the  computation  count 
for  TLM  analyses  increases  faster  than  the  fourth  power  cf  the  linear  mesh  density,  these 
accelerating  features  enhance  our  ability  to  model  complex  structures  to  a  much  greater 
extend  than  the  mere  memory  size  and  speed  of  the  computer.  Procedures  for  fine  tuning 
of  wall  positions  have  also  been  described. 

Fliture  time  domiun  CAD  systems  will  most  likely  employ  dedicated  parallel  proces* 
sors  configured  in  a  3D  arr^.  F\irthennore,  the  specific  nature  of  discrete  time  domain 
algorithms  requires  optimization  and  synthesis  procedures  different  from  those  employed 
in  traditional  frequency  doroeun  CAD  tools.  These  include  the  implementation  of  mov* 
ing  boundaries  for  geometrical  tuning  during  a  simulation  as  well  aa  numerical  synthesis 
through  reversal  of  the  TLM  process  in  time.  It  is  conceiv'able  that  at  the  present  rate  of 
progress  in  time  domain  modeling  these  procedures  will  equal  or  suipass  the  capabilities 
of  frequency  domain  CAD  tools  in  the  next  decade. 

REFERENCES 

[Ij  S.  Akbtarzad  and  P.B.  Johns,  ^Solution  of  MaxwelFs  Equations  in  Three  Space  Di* 
mensions  and  Time  by  the  T.L.M.  Method  of  Analysis,”  Proc.  lEE,  vol,  122,  no.  12, 
pp.  1344*1348,  Dec.  1975. 

|2J  K.S.  Yee,  “Numerical  Solution  of  Initial  Boundary  Value  Problems  involving  Maxwell's 
Equations  in  Isotropic  Media,”  IEEE  IVans.  Antennas  Propagation,  vol.  AP*14,  no. 
6,  pp.  302-307,  May  1966. 

[3]  P.B.  Johns,  “A  Symmetrical  Condensed  Node  for  the  TLM  Method,”  IEEE  TVans. 
Microwave  Theory  Tech.,  vol.  MTT-35,  no.  4,  pp.  370-377,  April  1987. 

[4]  Z.  Chen,  W.J.R.  Hoefer  and  M.  Ney,  “A  new  Finite-Difrcrence  Time-Domain  For¬ 
mulation  equivalent  to  the  TLM  Symmetrical  Condensed  Node,”  in  1991  IEEE  Inti. 
Microwave  Symp.  Dig.,  pp.  361*364,  Boston,  Mass.,  June  11-13,  1991. 


6 


[5]  J.  Nielsen,  W.J.R.  Hoefer,  “A  Complete  Dispersion  Analy^s  of  the  Condensed  Node 
TLM  Mesh,”  in  4th  Biennial  IEEE  Conference  on  Electromagnetic  Field  Computation 
Dig.,  Toronto,  Ont.,  Oct.  22-24, 1990. 

[6]  P.P.M.  So,  W  J.R.  Hoefer,  “3D-TLM  Time  Domain  Electromagnetic  Wave  Simulator 
for  Microwave  Circuit  Modehng,”  in  1991  IEEE  Inti.  Microwave  Symp.  Dig.,  pp. 
631-634,  Boston,  Mass.,  June  11-13, 1991. 

[7]  P.P.M.  So,  Eswarappa,  W.J.R.  Hoefer,  “A  Two-dimensional  TLM  Microwave  Field 
Simulator  using  New  Concepts  and  Procedures,”  IEEE  TVans.  Microwave  Theory 
Techniiiuos,  vol.  MTT-37,  no.  12,  pp.  18n-1884,  Dec.  1989. 

[8]  M.A.  Morgan,  Editor,  '‘Finite  Element  and  Finite  Difference  Methods  in  Electromag¬ 
netic  Scattering,”  PIER  2  Progress  in  l^ectromagnetics  Research,  Elsevier,  1990. 

[9]  C.E.  Tong,  Y.  Fidiuo,  “An  Efficient  Alg<mthm  for  Transmission  Line  Matrix  Analysis 
of  Electromagnetic  Problems  using  the  Symmetrical  Condensed  Node,”  IEEE  Ttans. 
Mi*' .'owavc  Theory  Techniques,  vol.  MTT-39,  no.  8,  pp.  1420-1424,  Auj..  1991. 

[10]  J.L.  Dubard,  D.  Pompei,  J.  Le  Roux,  A.  Papiernik,  “Characterization  of  Microstrip 
Antennas  using  the  TLM  Simulation  Associated  with  a  Prony-Pisarenko  Method,” 
Inti.  Journal  of  Numerical  Modelling,  vol  3,  no.  4,  pp.  269-285,  Dec.  1990. 

(U)  P.B.  Johns,  “IVansient  Analysis  of  Waveguides  with  Curved  Boundaries,”  Electronics 
Letters,  vol.  9,  no.  21, 18th  Oct.  1973. 

(12)  R.  Sorrentmo,  P.P.M.  So,  W.J.It  Hoefer,  “Numerical  Microwave  Synthesis  by  Inver¬ 
sion  of  the  TLM  Process”,  in  21st  European  Microwave  Conference  Dig.,  Stuttgart, 
Gennany,  9  •  12  Sept.  1991. 


7 


© 


INCIDENCE 


SCA'n'ERING 


•  1/2 

t 

.1/2-^.— *1/2* 

i 

1/2 

•  •  • 


(a)  (b) 


J 

♦- A/ 

,1V 

(e) 


‘-1/2V 

I/2V 

1/2V 

I/2V 

1 

(f) 
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Corresponijing  signo<-  Mow  (wave-flow)  represcntolion 
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For  deriving  ^  ^rom  continuous-domain  MD  circuit,  the 
corresponding  discrete-domain  HP  circuit, appi)' substitution' 
If  simple  sampling  is  cArried  out  in  originejl  coordinates^  t  • 
^  ^ (i  0, 0,Xi){ •  },  T --  spatial  shift 

^  A(0,±T,0,T^){.|;  T^=+ime  shift 
iA(0,0,tT,T,)^.i,  T,  =  T/v- 

1  A(o,o,o,aT,){  .  i 

Similarly  (^r  J>yi(£"-4')  • )  =  ((f'- 4)  •)  ■ 

2- More  efficient  (densest  hall  packing !),  but  less  eas^^ 
if  simple  sampling  is  carried  out  in  rotated  coord  mates,  t( 


Ps 


\  c 


Series  conneclion  of 
n  ports  (numbered  v=  I  to  n). 
Rv  =  port  resistance 

Vi.V;.-  •  .Vn=0 

*1  =  *2  =  ■  =  ’n 

Ov”  *  ^v*v  •  bv”'^v"^v*v 


Thus,  have  3n  equations  in  4n  variables 
Eliminate  all  v^,  ly  .  Solve  for  the  by  •- 

by  =  Qy  -  YyC  Qj  •  02*  *.0^),  Yv"  *  ^2*'  *  ^n) 
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Y^=2R/(R|.R2.  .  Rn),  Y,  *Y2  — •  •Yn=2 


Can  choose  one  port  as  dependent  port. 

e  g.  v=  1 :  Yi  =  2  -  Y2  ■Y3 "  '  '  "  Y^-  Thus,  need 

n-1  multipliers  (not  n^)  =  number  degrees  ol  freedom. 


Recent  Developments  of  in  Numerical  Integration 
of  Differential  Equations 
Wolfgang  Mathis 
UnivcrMty  of  Wuppertal 

Abstract 

Nixmerical  Integration  of  differential  equation  is  a  standaird  discipline  in  numer¬ 
ical  mathematics  and  basically  for  simulating  dynamical  systems  m  all  areas  of 
engineering.  In  dependence  of  the  kind  o^  modelling  dynamical  systems  \rill  be 
described  by  ordinary  or  partial  difierentiil  equations.  In  this  paper  we  restrict  us 
mainly  to  the  former  case.  The  most  general  type  of  ordinary  differential  equation 
has  the  implicit  form 

F(x,i,t)  =  0,  (1)' 

and  is  called  differtniiaUalgehraic  tquatton  (DAE).  By  means  of  a  suitable  trans¬ 
formation  the  explicit  dependence  of  t  can  be  dropped  If  F  is  solvable  for  x 
globally  we  obtain 

X  =  J(x);  (2) 

we  refer  this  type  as  ODE  Because  these  equations  possess  a  manifold  of  solutions 
rather  than  a  unique  solution  additional  conditions  must  be  prescribed  to  x.  If  x 
is  determined  at  one  time  to  this  situation  is  called  initial  valued  problem  (IVP) 
whereas  conditions  in  different  time  points  are  called  boundary  valued  problem 
(BVP). 

In  order  to  simulate  a  dynamical  system  we  start  with  a  system  description  of 
*yp«  (1)  (2)‘  In  dependence  of  our  interests  we  formulate  a  IVP  or  a  BVP.  For 

calculating  the  solution  x(f)  we  associate  a  convenient  difference  equation  to  (1)  or 
(2).  It  IS  obviously  to  replace  the  derivative  x  by  a  difference  approximation  with 
a  step-sire  h  and  to  construct  a  discrete  sequence  of  x(fn)>  The  quality  of  such 
approximations  will  be  characterized  by  consistency,  convergency  and  stability. 
The  theory  of  numerical  integration  of  ODE  (1)  and  the  art  of  its  implementation 
are  available.  In  this  paper  we  are  interested  mainly  but  not  only  in  multistep 
methods.  In  most  applications,  e  g.  mechanics  and  electrical  engineering,  most 
dynamical  systems  are  represented  by  DA£*s  in  a  natural  manner.  For  this  reason 
we  discuss  the  main  aspects  of  to  the  theory  and  implementations  of  numencal 
integration  methods  for  DAE’s  and  remark  that  the  essential  results  are  developed 
during  the  last  ten  years. 

Furthermore  wo  will  discuss  the  problem  of  step-size  control  and  switching  between 
different  integration  algorithms  (order  control)  &om  a  control  theoretical  point  of 
view;  this  part  includes  also  some  results  worked  out  m  our  corresponding  project. 
We  illustrate  this  mater.al  by  means  of  some  examples  from  circuit  simulation. 
The  hnal  section  contains  something  about  the  problem  of  rounding  errors.  This 
is  an  essential  subiect  because  the  development  of  algorithms  (operations^  and  the 
characterization  of  their  properties  will  be  discussed  in  the  real  numbers  JR  and 
m  other  sets,  e  g.  which  are  constructed  in  a  ‘vertical  manner*.  In 

classical  numencs  we  choose  a  suitable  finite  set  (e.g.  Coating  point  numbers  F) 
of  IR  and  associated  operations  and  construct  the  ‘higher*  sets  and  operations  in 
a  ‘vertical  manner’m,  too  Therefore  it  is  not  clear  that  the  numcncaf  algonikm 
(implemented  in  F)  works  in  the  manner  as  the  algorithm  in  SI,  To  circumvent 
this  problem  it  seems  to  be  useful  to  apply  a  well-defined  antbmetic  (Kuhsch- 
aritbmetic)  and  well-adapted  algorithms.  This  approach  is  close  related  to  the 
wave-digital  filter  method  of  Fettwds. 
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Caiie  step  site,  that  is,  A  0  (A  /•  0) 

•  Lieear  test  ODE  (l■dlBealleaa]  case) 

a?  (•» 

•  LMS  (*)  IS  absolot  stabll  at  AA  ,  it 

P(i)-(AAHO-0 
kaaroeis|r,)  <  1  I  •  J,  .k 
Denotation] 

$  a  {AA  C  C|  LMS  lorihod  is  aUjlst  itabd  fecAA  } 
it  called  Stability  Region 

•  Stability  region  of  the  test  ODE  t' m  >Ar, 
lolsitoa  ({A)  «  ('re  where  i  a  AA 

•  Stability  region  of  a  DE  ftom  the  test  ODEi 
approx  mIuCiob  ri  a  ^{i)ri  wc'ri 


l» 


W 


•  Ideal  Condition 

A.aC 

Example  Ctaah  Nicotson  Mtthod 


•  Daklquist's  A  Stability 

Aaalyitcal  drerrastag  awl«iioa>(ODC)  — •  decceasiag  soIvii«as(DC) 

S„3C‘ 

Example  laapliol  Euler  Method 


*  Unfortuoetely  A  SlaULly  trsiricii  Accuracy 
Tbeorem  (Daklquist) 

Espbcii  IMS  mtihodi  caaaol  be  A  ttalil 

[•pUil  LMS  wrlhodi  uilh  vcdrt  fiS  caaaot  be  A  slaUI 


b  saall 


si.l'ility  d«terfflio«s  fa  for  'SdfT  ODE** 


4  S  ‘Stiff*  ODB 


Exampl« 

/'  .  -XU  -  p(i))  -  ,(0)  ••  », 

-  S<J»tK)n  r(l)  m  (t,  -  +  p(0 

-  Ac(or»<v  PciaoftA]  LTE  f<»t  expliot  Ealti  s]*lk«ds 
f  <)  Uitul  Ttaiu«M 

b  imaU 

I  >  ]  Slow  vanatiMi  with  p 

-  Stabilty  Stadyiii  lb«  |l«bal  «rre?  <.  *  >((.)-{. 

f..,«|l+a)«  +  irE, 

t»  I*  unpbCtd  «»)«» 

-2<a<o 

Subdity  leilrKl*  b' 

•  Th(u  tim*  iMtitilt  ef  ODE'i  *rf  <>U<d  ‘ilifC 


AaODEuuid  to  be  cliff  oa  ii  ta«r«  euCc  a  eoaiposcst 
«f  a  coialMa  that  vanet  Ur{e  eonpaced  t»  1/7 
P«r  baeax  itoe  lavanaat  ODE  s  x'  ■  Ax 
Tb«  ODE’e  are  ‘ct.ff',  1/ 

»ax|A,(A)|  >  iiua|A.(A}| 


a  Remark  The  lact  deCaitioe  it  lavabd  hi  blear  time  vinanc 
ODE-c 
Example! 

+  IJ<ot*6/  +  9/Jci»I2l  \ 

\-12tiB*6l  +  9/Inal21  -1  - «cia’6< -««» J2t j  * 

Ei|raviltei  laid  10  (coatlaat') 

ExMt  Soliiiicn 


*{0  •  Cje**  ^ 


£616/4  21 
2c6i6t  -1 


>6/-2c6s$/\ 
in  6/  4  CM  6// 


latttpteUiion  ezp(>/]  aad  »p(-10<)  are  aot  laclsded 
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«  Th*  fiiilia]  iranclent  it  not  'ttlfT',  beecute 
|AA|  It  tmeU 

Example  Naa  d«f  fvl  F^gaii<n 


•  Detitn  Approaefaet  for  ‘nonitirr’  and  ‘ttlft’  ODE'e 

-  Clvb*l  Error  E^yaiioa 

/..I  •5,r.+  ITE, 

■  DetiKii  |oa]  la  noitiiff' racei 

IJ  C,  iniatl  at  potaiUe 

-  Detitn  goal  la  tliff  caert 

S,  imall  at  poitiUe 
Piiif  t7E.  lot  miaimal 


Diffvfenl  Stability  Concept!  for  Stiffnew 

ODD 


hmNImOki*  CftAleiii  eijppruiw 


ACtOlgU  SuC-le)! 


ISuMar  CaMncti«>)Cmi|iu 


! 


Further  Characteniations  of  ‘Stiffness’ 

•  L»rje  Step  Sue  h 
■  ImfJicit  aieik«>ie 

^v.ej  opg 

-  -*  B«..»e  of  !AA|  >  I  ..«p,o 

iecreuiel  ' 

(S»Bdl,erttSki<litr,»B,  l%j) 

e  SllffoeM  requires  Newton  type  Methods 

-  .-.ru„  PO.M  for  h.wto.  typ,  method. ..  seeded 
e  PrrdKior  neikod  deiermuei  st«rtis|  poiti 

.  Pr^«,or  Coneetor  d.ffer.ee.  prov^Je,  s  i.TE  e.t™. 


5  Remarks  to  the  Implementation  of  ODE-Solvers 


Problems  for  ImplementlB*  a  IMS  family 
e  Soluble  P^ysomis]  Repreteotation 
•  Etficieat  Step  Site  Coairo! 
e  EflkieAt  Order  Cotirol 

-  Very  Usefa]  Ceaerpii  from  Coattol  Theory 
A»  Eiamplc)  Rja(  Modolsior 


UD> 


Furih«  Rem»rk»  to  Impl«m«nt*tion» 


‘  *  Suilsm  I)«tfr(>o» 

I  .  EfGotii  NcBltM,  S„|.„  lyp,  Mttkiodi) 


f 


i 
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6  nevwed ‘Desc.iption  Equation  for 


•  E*:J«rtODEi»f«#etMluf»Uy  fof  CiKuil  Sim«l*tieB 

'  S«  EtBBiplr  ,s  SfttK.«  I  impUct  ODE  d.rrctly  ob[*,fl,d 

•  Equaiioet  (D^E) 

-  l®t>l»otODE(«DAE) 

F(K.x'.0-e 

-  StfBi  fxpjifrt  DAE 

x'-flx.y) 

o»f(*.y) 

•  DcctrKBl  Nttwofli  trt  deirribcd  by  t  Dutsrr  of 
•  Al|rbr«tc  oqsaueni 

(knAWll..  i,,., 

»«>•*.  unttt  tk*r*ti*rit»t)OB) 

~  DifferfBiiU 

(!.•.♦«  Btd  ooobtMr  f.pMitcr  *«d  ladiKtot  cbt»ctrr,i« 


DwctiUin  N«t*efki  by  Srn.i  DAE'* 
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fotaOMitt* 
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i.  Numerical  Solutions  of  DAE’s  9-  Summary  and  Outlook 
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Definlllon  of  CA 


Cellular  Automata: 
Applications  and  implementation 


Lolhar  Thidc 

I>*hrs(uh]  fur  MikrotUltronik 
Unntrsitll  das  Sairlandti 
D  (600  Saarbrucktn 


Formal  Delinitton 


{\\I  /€l  deP))) 

1  site  index 

1  index  domain  of  CA 

(;  aibiirary  function  0  5'^*  -♦  S 

P  neighborhood,  e  g  T  =  J  6  Z”  A  [I  /|!-^  <  >■} 

S'  sc!  of  states,  i  c  a(7,l)€5 

'  set  of  configurations,  i  c  E  =  S'* 

‘I’  global  mapping,  i  e  4’  E  E 

set  ol  configurations  generated  after  i  iterated 
applications  ol  i  e  G  E, 


I 

Discrete  m  space:  CA  consist  oi  a  discrete  lattice  of  sues 
Discrete  iit  time*  CA  evolve  m  discrete  time  steps 

Discrete  states*  Each  state  takes  on  a  finite  set  of  possible 
values  'i 

Homogeneous  All  cells  are  identical  and  are  arranged  in 
a  regular  way 

Synchronous:  All  cel!  values  are  updated  in  synchrony 

Deterministic*  Each  cell  is  uixlaled  according  to  a  fixed 
determmisiic  rule 

Spaiiatl)  local*  The  nile  depends  only  on  the  values  of  a 
local  neighborhood 

Temporally  local.  The  rule  depends  only  on  values  foi  a 
fixed  number  of  preceding  steps 


_ Related  Models _ 

Partial  difTcreniial  equations;  space,  tune  and  sue  values 
are  continuous 

Finite  difference  equations:  sue  values  are  contuiuous 

Particle  models:  panicles  have  continuous  positions  and 
velocities 

Neural  nelnork  models,  connection  patterns  aie  aibitiaiy, 
Site  values  are  continuous,  updates  are  asyncluonous 

Cellular  neural  networks*  sue  values  are  continuous 

Iterative  arrays:  different  purpose  than  CA 

Array  processors:  sites  can  store  extensive  infonnaiion 


Applications  of  CA 


Physical  Modeling 


(  ipdtation  Theory 

•  self  rcproduvijon  (J  v  Neumann  1949.  Conway  1970) 

•  fonnal  language  (S  Wolfram  1984) 

•  cljssificanon  (S  Wolfram  1984/1985) 

Biological  Modeling 

•  self  reproduciion  (J  v  Ncumarm  1949) 

•  evolution  (S  Ulam  1918) 

•  mcdeli  of  memory  (M  Minsk)  ’969) 

Physical  Modeling 

•  calculating  spaces  (K  Zuse  1940) 

•  hydrodynamics  (J  Hardi  1973,  U  Trisch  1986) 

•  growth  mechaiusms  (J  D  Gunton  1983) 

atiem  recognition  (K  Preston  1984) 

•  spin  models  (M  Cieutz  1980) 

•  wave  models  (H  Chen  19S8) 


JB  Salem,  S  Wolfram  'niermodymamics  and  Hydrodyna* 
niics  with  Cellular  Automata  In  llicory  and  Applications 
Cullulai  Autuifiula  WuilJ  Svieiitiriv,  1987 


Motivation  1 

,•  phenomena  are  often  desenbed  by  (nonlinear)  (paniai) 
differential  equations,  two  different  approaches  to  s.i.n<- 
laiion. 

1  •  diffvicnce  equation  on  a  macroscopic  scale 

•  discretization  ui  time  and  space 
-  •  microscopic,  discretized  mode’ 

•  large  number  of  similar  components  with  local  m 
tcraction* 

•  functional  homogeneity  reflects  space  and  time  invari¬ 
ance,  locality  reflects  finite  speed  of  information 

•  cellular  automata  arc  simple  to  program  and  amenable  (o 
parallel  processing 

•  studies  of  collective  phenomena  possible  (turbulence, 
chaos,  fractality,  ) 


_ Fluid  Dynamics _ 

(Hardy  Pawn  Pomeau  1972.  Fnich  Hulacher  Pomeau  1986) 

fluid  incompressible,  absence  of  external  forces  Navier- 
Stokes  I^uation  (NSE) 

^  +  VAV  =  --Vp  +  lV*l' 

Ot  p 

V  velocity 

p  (constant)  density 

V  viscosity 

fictitious  microscopic  model' 

•  as  simple  as  possible  dynamics  (not  necessarily  follo- 
wmg  Hamiltonian  equations  for  mteracting  panicles) 

•  reproduces  NSE  on  r  croscopic  level 

particles  and  their  scattcruig  ore  moulded  by  reversible 
CA  rules 


B  Salem.  S  Wolfram  T!.emiodynamics  and  Hydrodyna¬ 
mics  w.ih  Cellular  Aulomaia  In  Theory  and  Appl, canons 
of  Cellular  Amomala  World  Sciemific,  1987 


- - - - 


;U 


Flow  pisl  an  obstacle  (from  Salem  and  Wolfram) 


- - - Laltiee  Gas _  ■ 

•  Panicles  move  on  a  lattice  and  satisfy  certain  symmeliy 
requirements  Movmg  and  scattering  by  reversible  rules 

•  Dcnvation  of  macroscopic  behavior 

•  molecular  level  motion  is  reversible 

•  Imelic  level  nonequilibrium  statistical  mechanics 

•  macroscopic  level  conlmuum  approximation 

‘  Wave  equations- 

•  foi  small  penurbaiions  from  equilibrium-  linear  elastic 
propenics  of  lattice  gas 

•  piopagaiion  of  a  distuibatioii  is  governed  by  wave 
equation 

•  study  of  wave  interference,  reficciion.  diffraction,  re- 


Modeling  of  Wave  Equations 

(H  0>cn  1988)  ‘ 

•  Linear  wave  propagation 

at- 

•  Invariants 

energy  ~  J  {{ji) 

momentum  P  =  '^J  {  (^)  ^“  "  (It) 

•  Concept 

•  two  kinds  of  photons  propagating  on  a  lattice  <r  =  f . 

.  =  -?  I 

‘  number  of  photons  with  quantum  a  ata 

panicular  sue  x  and  time  i  moving  wuh  velocity  Ca  ' 


•  Hyugens  pnnciple  any  spatial  point  can  be  drought  of 

as  a  new  wave  source  with  intensity  | 

s=  — u(a:,t)  ' 

m  I 

•  decay  rate  of  source 

S(l,  (  +  1)  =  S(jr,  0+ 5!^  <r  (At" (i,  ()  -  At^(r  +  Ca,  0) 

e,c 

+  1)  s=  u(x,t)  -  g{x,t  +  1) 

•  continuous  Imcar  wave  equation  is  recovered  after  ma* 
king  an  ensemble  averagmg 

•  result  can  be  converted  to  a  finite  difference  eouaiion 


H  Clicn,  S  Chen.  G  Doolen,  Y  C  Lee  Simple  Lat 
iiie  Gas  Models  for  Waves  Complcv  Systems  2  (1988) 
259-267 


Implementation  of  Cellular  Systems 

(SY  Kgnj,  P  DtwiWe,  E  Dcprnitie  S  M«rk«) 

Specificalion 

•  Cellular  automata 

(M  /el  fl(/,0-d(M/-d.t-l)  deV})) 

•  Cellular  Systems 

(11/  /€l 

ai(/)‘=<^l({ni(/-d}  dePii},  , 

{ar(/-£i)  d€Dri})|l 


) 


av(/)  =  d'i({oi(/-d)  d€l>ir}.  . 

dePnl) 


inUrfcrcuct  I'tlUfA  o'  iVf  v>«v(  atni  Iilvic  »  a  two- 
il  mcKi>or&'  h*vc  U'.lKc  (uCoiUc  v'll  cki'crimcrt  at  a  (itcfl  ruatar; 


Program 


(II  0<i<2a<>0 

l,i  -  -  D)  II 

Q2(^0  =  -  1)) 

) 

Juced  dependence  graph 


Target  Archileclure  j 

•  dcdjcaied  hardware  (CAM-6.  CAM-8  (Toffoh))  | 

•  coarse  gram  parallel  systems  (MIMD,  Transputer)  ' 

•  fine  gram  parallel  systems  (SIMD,  Coiuicciion  Machini 

Mapping  criteria  i 

•  communication  «— <  computation 

•  consideration  of  pipclmed  antlimetic  units 

•  consideration  of  finite  resources 

•  suited  for  automatic  compilation 


Applications 

•  neural  networks 

•  Iterative  arrays 

•  cellular  automata 

•  solution  of  PDC 

•  systolic  anays 


logeihw 


Implementations 


Mapping  Problems 


n%uralion: 


Partitioning  (limited  resources) 


Scheduling  (pipelined  arithmetic  units) 


0(i»r»t«ns 


Specification  of  embedding  (soflMarc) 


Partitioning 


Partitioning 


Projection. 

•  reduce  dinieiuion  of  DO  by  1 

•  affine  projection  of  DO 

ultiprojecdon' 

•  reduce  diincnsion  of  DG  arbiuanly 

•  loop  control  or  flow  control 

•  inai^h  I/O  rate  and  computation  rate 

DG  Sequencing  Signil  flow  graph 

=  rOrS-g:: 


Actne  Clustering. 

•  match  given  number  of  cells 

•'  match  gi\en  dimension  of  array 

•  combine  muUiprojcction  and  passive  clustering 

1  ocal  Sequential  /  Globa)  Parallel  (LSGP)i 


Processor  array  Clustered  anay 


Passive  Clustering' 

•  make  inefficient  array  efficient 

•  nuke  use  of  pipelined  units 


Local  Parallel  /  Global  Sequential  (LPGS)* 


Processor  array 

o®4-6-e-B 


Clustered  array 


Hierarchical  Transformations 

ling  of  Iteration  space  using  tiling  matnx  P 


Hierarchical  Representation  of  Algorithms 

•  Each  node  represents  a  program  containing  fimctions  delinei 
in  lower  hierarchical  levels  | 


•  iniroduLiion  of  cennmals 


program  iransformaiion 
always 
{11/  /€l 

ahva\s 

(1]J,  A,*r  JiJAh’^KA'r^V 
zlJ",  A']  -  /{z[7  Fi  iw  F-)]} 
if  J-d- 


•  There  are  two  basic  program  schemas 

a  Straight  line  code  (set  of  assignment  siaiemcnis.  nonregu’ai 
dependence  graph) 


b  Regular  algorithm  (regular  dependence  graph) 


Hierarchical  Compilation  Strategy _ 

seiiucnce  of  provable  correct  program  transformations 
specification  of  paramcieriwble  basic  iransformaiions 
s  ptimization 


,  spscification 

.  mr 

!  program 


_ Mathematical  Tools 

•  opeiations  on  index  sets 

•  tiling  of  Iteration  spaces 
♦.  mieger  linear  algebra 

•  linear  and  integer  linear  programming 

•  geometry  of  numbers 

•  combinatorial  algorithms  on  periodic  graphs 


Compiler  For  Cellular  Systems 


Discrete  Time  Domain  Modelling  of  Electromagnetic  Fields  and  Networks 
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ANALYSIS  OF  NONLINEAR  MICROWAVE  dRCUITS  VIA  THE 
TIME-DOMAIN  VOLTAGE-UPDATE  METHOD 

H.D.  Foltz,  ‘JH.  Davis,  and  T.  Itoht 

Although  direct  transient  time-domain  solutions  of  circuit  equations  (for  example,  SPICE-type 
programs)  are  sometimes  employed  to  analyze  microwave  systems  containing  nonlmear  devices, 
steady-state  methods  (for  example,  harmonic  balance  techniques)  are  preferable  in  cases  involving 
bgh-Q  resonant  cucuits  and  other  narrow-band  structures,  fm  which  steady-state  methods  have  a 
significant  advantage  in  computation  tune. 

Many  steady-state  nonlinear  techniques,  such  as  (1)  the  common  minimizing-or-optimizing-bascd 
piecewise  haimonic-balance  method  and  (2)  the  relaxation-based  voltage-update  method,  operate 
piimarily  in  the  frequency  domain.  Unfortunately,  most  nonlinear  solid-state  devices  are  most 
easily  treated  in  the  time-domain.  Analagous  steady-state  techniques,  based  on  discrete  time 
samples,  can  be  foitnulated;  (1)  the  waveform-balance  technique,  which  is  related  to  the  piecewise 
haimonic-balance  method,  and  (2)  the  time-domain  voltage-u^ate  technique,  which  is  relaxation- 
based.  In  this  paper,  the  latter  technique  will  be  examined  in  more  detail  and  compared  to  more 
conventional  approaches. 

In  the  time-domain  version  of  the  voltage-update,  time  samples  representing  the  steady-state 
voltage  waveform  are  applied  to  the  nonlinear  device(s).  The  resulnng  current  samples  are  applied 
to  the  linear  pomons  of  the  system,  leading  to  a  new  set  of  voltage  Samples.  Relaxanon  parameters 
are  then  applied  to  determine  the  starting  samples  for  the  next  iteranon,  and  the  process  is  repeated 
until  convergence  is  obtained 

Voltage-update  techniques  have  a  marked  advantage  over  most  other  approaches  in  simplicity  and 
speed  per  iteration,  when  applied  to  problems  in  which  the  frequency  is  known,  such  as 
amphfiers,  frequency  multipliers,  and  tmxers  They  can  also  be  applied,  with  some  modificanons, 
to  vanable-frequency  problems  such  as  oscillators. 

Strategies  for  extending  the  range  and  speed  of  convergence  for  the  relaxanon  procedure  will  be 
discussed,  along  with  the  relanonship  between  frequency-domain  and  time-domain  relaxation 
parameters.  The  results  of  several  representative  appUcations  mvolving  negative  resistance  devices 
and  SIS  junctions  will  be  presented. 


•  The  Umversity  of  Texas.  Eiectncal  Engineering  Research  Laboratoiy,  Austin.  TX  787 12 
f  ^versity  of  Cahfomia,  Los  Angeles,  Department  of  Eiectncal  Engineenng,  Los  Angeles,  CA 
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PERIODIC  POWER  COMBINER 


•  D*vie*i  oomintljy  t»U*v*vele8tlh  aptrt 

•  Load  ifflpodascc  >  device  Impedaece  /  N 


•  Oevicet  resonated  van  indueUve  siobs  ti 
be  resistive  at  oselltation  fredueney 
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LOCKING  BANDWIDTH 
VERSUS 

INJECTION  LEVEL 


S  .f  I 

^  4.0  OZ  M  If  AS  II 

IrOiaiON  CURRENT  (BA) 


Multiple  diode  combiners  require  higher 
injection  currents  to  obtain  a  given 
locking  bandwidth 


INJECTION-LOCKED  POWER 
COMBINER 


•  model  INJECTION  WITH  CURRENT  SOURCE 
AT  END  OF  ARRAY 


■  ''O'*  frequency  UPDATE 
FREQUENCY  IS  KNOWN 

■  USE  TIME  DOMAIN  UPDATING 

■  yi^F^JES  OF  CANCELLING  RESISTORS 
*X  device  port  to  improve  • 
CONVERGENCE 


OSCILLATION  FREQUENCY 
WITH  RESONANCE  FREQUENCY  OF 
INDIVIDUAL  DEVICES 


»<  >(  U  111  IIZ  IM  III 
IN&tVOiUAL  DEVICE  *E$ONA^CE  (CHi) 


As  the  number  of  devices  In  the  combiner 
Increases*  the  oscillation  frequency  depends 
more  on  the  devices  and  less  on  the  periodidty 


solution 


VARIATION  OFTOTAL  OUTPUT  POWER 
WITH  RESONANCE  FREQUENCY  OF 
INDIVIDUAL  DEVICES 


variation  of  POWER  DICTRICUriON 
AMONG  DEVICES 


™™vALD.V.CtR«ON«.CE<0„.,  As .h. b  ch.„,.d  from .h. 

ctnter  frequency  where  the  structure  fs 
half-wave  periodiCt  the  distribution  of  power 
becomes  unequal 


EFFECT  OF  DEVICE  NON.UNIFORMITY; 

FREQUENCY  AND  POWER  OUTPUT 
VERSUS  DEVICE  AREA 

(THREE  DIODE  COMBINER) 


«j  S.0  14  xe  u 


AXCA  or  rERTVRBED  DEVICE  REUTIVE  TO 
OTHER  DEVICES  IN  COMBINER 


CONCLUSIONS 


Time  Domain  Voltage  Update  Method: 

-  Useful  alternative  to  other 
steady-state  methods 

-  Has  been  applied  to  wide  variety 
of  circuits 

-  Should  combine  easily  with 
electromagnetic  structures  analyzed 
by  time-domain  techniques 
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The  state,  output  and  non-linear 
equations  of  the  whole  network 
consisting  of  a  number  of  multiports 
is  written  in  the  form, 
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ESdent  Analytteal-Numerical  Xlo^Uag  0{t]ItTa>'W!de1)«a(i  Pnlied 
Plaae  Wave  Scattering  From  a  tugs  Strip  Offtiog 


t  Lawrence  Carln  and  Leopold  B.  Felten 

Weber  Research  Institute/  Electrical  Engineering  Department 
Polytechnic  Univenity 

1  Familngdale,NYH735 


Summary:  Ultra-wideband  (UWB)  pulsed  plane  wave  scattering  from  a  large  but  Mte  atrip 
grating  in  free  space  is  analyied  in  the  frequency  domain  via  decomposition  into  plane  wave 
spectra,  implemented  numerically  by  the  method  of  momenta,  and  then  inverted  into  the 
lime  domain  (TD).  To  make  this  procedure  practical  under  UWB  conditlonf,  doted  form 
■  expressions  are  derived  for  Interaction  integrals  involving  widely  separated  and 

testing  functions.  These  closed  forms  are  based  on  ajudicious  choice  of  the  basis  funetioaa, 

and  on  asymptotic  methods  for  reducing  the  integrals.  Although  large  separation  distances 
are  assumed,  the  eapressions  have  been  found  to  be  accurate  for  leparaisont  at  small  as  0.1 
wavelengtha  The  TD  self  terms  can  also  be  evaluated  efBdenUy.  To  test  the  frequenty 
•  domain  algorithm,  comparisons  are  made  with  available  deta  in  the  literature  for  ,  i 

currents  and  far  Geld  scattering  from  multiple  strips.  Newshort  pulse  TD  results  are  shown  '  1  I 

as  well.  ' 


Transient^^urr*nts^and;^fieW^^o^»ire^ntenna^jitl^iode^ 

N.  Scheffer,,  Telefunken  Systemtechnik  VR3  E51,  Sedanstr.  10,  7900  Ulm 
l^^JJuingricaljnethod 

The  electric  field  integral  equation  for  a  ware  antenna  can  be  written  in  matrix 
notation  for  the>time  step  to  as 

(2+2f).j,-E;+E;.  (1) 

The  unknown  quantitiy  in  <1)  is  the  current  distribution  !,■ 
of  the  antenna.  The  antenna  is  devided  into  N  segments. 

The  impedance  matrix  Z  is  time .independent. 

r  ft  °  ^ 

Zf  •  I  I  represents  a  diagonal  matrix  containing 

I  0  R»} 

the  nonlinear  loads. 

The  components  of  E;  are  the  tangential  components  of  the  incident  electric 
field  and  £f  are  the  tangential  components  of  the  scattered  electric  field. 
Introducing  the  admittance  matrix  Ya2*'  *nd 

(2) 

equation  (1)  can  be  written  as 

(E  +  Y-Zf)-!..:;.  (3) 


E  repesentr  the  unity  matrix. 

For  the  special  case,  that  the  antenna  is  loaded  at  exactly  two  segments  m  and  n 
with  and  k,  .  the  currents  in  the  loaded  segments  are 


- 

(y*.  A + jxy— + j)  - 


The  currents  in  the  unloaded  segments  ieim.n  are 


Y^IUK,  -  .  ( 6  > 


2.  the  tine  dependent  current  distribution 


The  above  method  of  calculation  is  used  to  examine  a  dipole.  The  dipole  contains 
two  diodes  and  is  centerfed  by  a  source  voltage  with  Gausssian  pulse  shape  as 
shown  in  Fig.  1. 

The  two  diodes  are  described  by  resistors  ... 

B  I  Oh  .  a  Sv 
"*‘*1  1004ft  <0 
B  /Oft  ,/;>o 

1004ft  .;?<o 

Fig.  2  shows  the  two  current  pulses  arriving  at  the  antenna  ends  after  passing 
the  diodes  almost  undistcrtedly.  The  arrows  mark  the  positions  of  the  diodes.  In 
Fig.  3  the  current  pulses  have  arrived  at  the  reversed  biased  diodes  after  they 
have  been  reflected  at  the  antenna  ends.  Fig.  4  shows  the  reflection  at  the  re¬ 
versed  biased  diodes. 


S^^^^Th^^ljec^jjCjiearfa^ 

It  IS  possible  to  reduce  the  fieldline  equation  irxSmi}  to  a  potential  .function 


:-f  I  # ■  V*’ f S'*  because  of  the  cylinder  symaetry. 

An  electric  fieidline  is  then  described  by  V-eowj  ,  is  derived  in  the  ap¬ 
pendix.  Fig.  5  shows  the  electric  fieldlines  located  on  concentric  spheres 
around  the  generator.  The  pulse  has  just  arrived  at  the  antenna  end.  In  Fig.  6 
the  reflection  at  the  antenna  end  is  considered:  the  fieldlines  are  now  located 
on  concentric  spheres  around  the  source  and  the  antenna  end.  Fig.  7  shows  the 
reflection  at  the  diode:  a  new  sphere  around  the  diode  is  visible.  In  Fig.  8 
multiple  reflections  for  a  late  time  step  can  be  observed. 


Calculating  Frequency  Domain  Data  by 
Time  Domain  Methods 

ty  M.  Drlilisr,  M.  Dohlus,  T.  Weiln  id 

Ttchnueke  UoehscLulc  DarMStadt  Tlicorie  £id(troma|;n^t»cil«r  folder  (FB  U) 

Aba^raet 

W«  sKow  th«  eal«uUtiea  iA  far  f!dd  pattarna  aad  Katt«iiii|  paramtim  by  aiaaa*  «( time  do* 
aaiA  method*.  la  order  to  o  jtaifi  a  inod«exciiaUoa»  the  c5iten«^tf(Ioa  of  the  diaercie  waveguide 
ci|eBvaI«e  problem  la  com!  mation  with  aa  adequate  procewag  at  the  boundarie»  it  ttaed. 

Far  FmU  IVaueforms 

lu  general  the  electric  field  c.  radiating  atzuciurea  can  l>c  written  for  largt*  dittancoa  a> 

(1) 

where  Jt  a  w^;u  u  the  wavr  mimber  and  r,6,d  denote  apherlca)  coordinates. 

For  the  calcnlatlou  of  the  far  held  transform  one  nroda  to  (ictcrmlne  the  eomplec  time 

harmonic  fuld  Amplitudes.  can  be  detcrmiuvd  hy  a  convolution  of  Greens  function  and 

tangeiitiAl  electric  and  magnetic  amplitudes  on  a  closed  surface  surrounding  the  radiating  structure. 

?(«,«)  =  jf/"''' {''*(&  x(i!x£))-  h-  x(>Tx£)|iM  ,  (2) 

The  fctnpuiaiion  of  the  complex  field  Amplitudes  by  a  Fast  Fontler  ’lYAniforin  ff(|uites  aampling 
and  storage  of  th«  tangential  electric  and  magnette  field  valu««  at  the  integration  surface  and  U 
only  feasiMe  for  small  mesh  »im  Therefore  these  ampUludee  are  obt.\lned  by  using  time  harmonic 
excitation  sources  and  a  direct  sampitag  of  the  harmonic  fields.  Another  way  is  to  perform  a 
monochromatic,  aingle  frequeney  fouiter  transform  of  the  lime  domain  fields. 

Filter  Dwslgn 

'fh«  transvcrcal  eliciric  held  iu  a  waveguide  ts  desortbed  for  the  time  harmonic  ca^e  by  the  mode 
expansion 

which  can  be  formulated  for  general  time  dependency  with  double  coiivolutiona; 

¥ 

With  /VO.r)  being  the  Inverse  Fourier  transform  of  exp(7,.(/w).  The  knowledge  of  4»(;w)  = 
&0>^)  ^  I*  esveritlal  for  the  calculation  of  the  tangential  boundary  field  £{n,2r,0tybr)  at  the 
interface  r  s  0.  The  stimulation  w.Ave  is  kstown,  but  ha^  to  be  calculated  from 

£.(jw)  s  ■  .  (^) 

£«0<*’)  be  obtained  by  mode  expaiision  of  the.  Held  ui  the  plane  t  s  St.  Tlie  modes 
and  the  projiagation  eoiuiunts  are  calculatvd  by  a  2D  eigenvalue  «o)ver.  The  case  of  frci^ueney 
dependent  ^  U  problematic  becaoK  the  expausion  of  sampled  fields  U  only  possible 

in  the  steady  state.  In  the  other  ease  is  yielded  directly  and 

A,(()  =  11-  +  BAt)  (6) 

can  be  calculated  by  single  convolutions.  This  method  is  applicable  for  hotnogeaeoualy  filled  wa\eg> 
uides  and  for  modes  with  wesk  frequency  dependency.  A  simpUneatiou  of  the  algorithm  is  possible, 
if  the  convolution  can  be  .approximated  in  the  desired  ficqueney  range  by  a  low  order  digital  HUer. 


TIME  DOMAIN  ANALYSIS  OF  INHOMOGENEOUSLY  LOADED 
STRUCTURES  USING  EIGENFUNCTION  EXPANSION 


MlCBAl  MR020WSKI 

Technical  University  or  Gdansk,  Telecommunication  Institute 
80-952  Gdansk,  Poland,  tel.;  +(048  58)  472  549,  fax:  +(048  58)  47  lo  7l 


At  present  there  are  two  algorithms,  namely  the  TLM  and  FD>TD,  which  are  used  to  solve 
Maxwell’s  equations  in  time  domain.  In  this  contribution  we  shall  present  new  methods  which 
may  broaden  the  range  of  options  available  in  time  domain  analysis  of  2«D  and  3*0  structures. 
A  wave  is  treated  as  a  superposition  of  tigenmodes  (eigenfunctions)  of  the  homogeneous  Laplace 
equation.  An  inhomogeneity  in  the  structure  perturbs  the  field  and  causes  the  coupling  of  eigen* 
modes.  Eigenmodes  are  chosen  so  that  they  fulfil  the  Helmholts  equation  either  on  the  entire 
homogeneous  domain  or  on  homogeneous  subdomains.  An  advantage  of  this  approach  is  that  it 
allows  to  obtain  time  domain  algorithms  which,  in  contrast  to  TIM  and  FD*TD  methods,  do  not 
exhibit  the  numerical  dispersion. 

Outline  of  time  domain  eigenfunction  expansion  algorithms. 

Based  on  the  concept  briefiy  described  above,  a  number  of  algorithms  can  be  proposed.  We 
shall  start  with  an  algorithm  called  a  complete  ^genfunction  expansion  (CEE),  let  us  consider  a 
set  of  coupled  differential  equations  reflecting  the  form  of  MaxweU  equations 

=  (1) 

where  £|,  £,  ue  Uneu  operators. 

In  the  FD-TD  elgoiithm  the  ebove  eqtietione  are  iliscretired  both  In  time  and  space.  In  the  CEE 
algorithm  the  discietiratlon  Is  onl"  In  time.  As  a  result  we  get 

/"  =  +  AtC,r  (2) 

The  unknown  functions  f,g  are  now  expanded  into  series  of  complete  set  of  orthonormal 
functions. 

/  =  ^«./.  ,  =  ^4,p.  (3) 

Expansion  function  are  defined  on  the  entire  domain.  A  sensible  choice  for  the  electromagnetic 
fields  is  are  the  eigenfunctions  of  Laplace  equation  with  suitable  boundary  conditions.  Substituting 
(3)  into  (2)  we  get 

Taking  the  inner  product  with  the  expansion  functions  results  in 

<1^  =  +  At  <  >  (5) 

=  &:■!/*  + At  <£3/",s.> 
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The  above  equations  can  be  cast  into  the  following  matrix  form 

a"=a”‘'  +  A<4i“*''’  (6) 

4"+>/»  =  i“-'/2  +  A)£a” 

where  s.  the  vectors  containing  expansion  coefficients  and  A  R  dense  matrices 

with  elements 

A,j=<C\gi,fi>  -fiiV  =<  >  (7) 

Another  version  of  the  eigenfunction  expansion  algorithm  is  obtained  if  the  discretization  is  in 
time  and  one  spatial  coordinate  and  the  expansion  Is  done  with  respect  to  two  remaining  spatial 
coordinate.  This  algorithm  we  shall  call  partial  tigenfunction  expansion  (PEE).  In  this  technique 
the  space  is  sliced  into  subdomains  and  the  fields  are  expanded  on  each  subdomain  (slice)  into 
series  of  local  expansion  functions.  In  the  PEE  method  one  obtains  a  set  of  equations  similar  to 
(6)  except  that  matrices  ^  R  sparse. 

Compared  with  the  FD*TD  method  the  CEE  and  PEE  algorithms  show  the  time  evolution  of 
the  expansion  coefficients  rather  then  held  components  at  nodes.  Such  an  approach  allows  one  to 
investigate  propagation  of  particular  modes  and  their  mutual  interactions.  Moreover,  in  contrast 
to  FD'TD  and  TLM  techniques,  both  algorithms  proposed  in  this  contribution  do  not  exhibit 
numerical  dispersion. 

EfHeient  numerical  implementation  of  eigenfunction  expansion  algorithms:  C£E*FFT 
and  PEE-FFT. 

One  drawback  of  the  CEE  and  PEE  algorithms  that  they  may  lead  to  higher  numerical  cost 
then  FD<TD  and  TLM.  The  CEE  involves  matnx  multiplication  hence-  assuming  that  expansion  is 
done  using  L  eigenfuncions,  the  cost  of  one  time  step  Is  of  order  0{L^)<  For  the  PEE  this  cost  is  lower 
as  the  matrices  involved  are  sparse.  In  the  FD*TO  and  TLM  method  with  nodes,  the  numerical 
cost  is  of  order  O(^).  Consequently,  eigenfunction  expansion  techniques  may  be  regarded  ac  an 
alternative  to  classical  time  dom^o  algorithms  only  when  ^  N.  This  condition  will  be  fulfilled 
in  slight^  and  moderately  perturbed  homogeneous  structures.  Nevertheless,  much  more  efhcient 
version  of  CEE  and  PEE  may  be  obtained  if  the  expansion  functions  are  sine  and  cosines.  Equations 
(6)  imply  that  at  each  step  one  evaluates  the  inner  products  <  >  and  <  >. 

For  sine  and  cosine  functions  the  inner  product  cvi  be  computed  m  a  very  efficient  way  using 
the  technique  described  in  (1).  In  this  technique  the  inner  products  are  computed  in  a  sequence 
of  inverse  and  forward  FFTs.  The  numerical  cost  of  such  computations  is  low  and  therefore  the 
overall  performance  of  the  CEE^FFt  and  PEE-FFT  algorithms  is  better  than  original  CEE  and 
PEE  methods. 

Conclusions.  New  algorithms  of  the  time  domain  analysis  of  inbomogeneously  loaded  microwave 
structures  have  been  described.  The  methods  proposed  are  based  on  the  expansion  of  fields  into 
complete  series  of  orthogonal  eigenfunctions.  The  resulting  equations  show  the  time  evolution  of  the 
expansion  coefficients  and  consequently  allow  one  to  investigate  propagation  of  separate  modes  and 
their  mutual  interactions  The  algorithms  proposed  in  this  contribution  do  not  exhibit  numerical 
dispersion  and  allow  coarser  time  discretization  than  the  equivalent  FD-TD  or  TLM  program. 

(1)  M.  Mrozowskt,  *'IE£M  FFT  •  A  fast  and  efHcient  tool  for  rigorous  computations  of  propa¬ 
gation  constants  and  field  distributions  lo  dielectric  guides  with  arbitrary  cross-section  and 
permittivity  profiles",  IEEE  Tnns.  Microvocvt  Thtorg  Ttch.,  vol.  MTT-39,  Feb.  1991. 
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The  Hilbert  Space  Formulation  of  the  TLM 
Method 

Peter  Russer,  Michael  KrumphoU  * 

Abstract 

The  Hilbert  space  represeatation  of  the  TLM  method  for  time  do* 
mein  computation  of  electromaipietic  Helds  and  the  algebraic  compu* 
tation  of  the  discrete  Green’s  function  are  investigated.  The  complete 
Held  state  is  represented  by  a  Hilbert  space  vector.  The  space  and  time 
evolution  of  the  Held  state  vector  is  governed  by  operator  equations  la 
Hilbert  space.  The  discrete  Green’s  functions  may  be  represented  by  a 
Neumann  series  in  spaco-  and  time-shift  operators.  The  Hilbert  space 
representation  allows  the  description  of  the  geometric  structures  by  pro* 
jectlon  operators,  too.  The  system  of  dlHerence  equations  governing  the 
time  evolution  of  the  electromagnetic  Held  In  configuration  space  is  de¬ 
rived  from  the  operator  equation  for  the  field  state  vector  In  the  Hilbert 
space. 


1  Introduction 

The  TLM  (tiwsmission  line  matrix)  metb(.d  developed  and  first  published  in 
1971  by  Johns  and  Beurle  is  a  discrete  time  domain  method  for  electromag¬ 
netic  field  computation  |1>2,3).  In  this  paper,  the  Hilbert  space  representation 
of  the  TLM  method  is  presented  and  applied  to  the  algebraic  computation  of 
discrete  Green’s  functions.  The  Hilbert  space  representation  is  a  very  general 
and  powerful  concept  in  field  theory  (4).  Whereas  in  the  electromagnetic  the¬ 
ory  Hilbert  s  »ce  methods  ace  mainly  used  for  solving  the  field  equations  as 
for  example  in  the  moment  method  (5],  in  quantum  theory,  the  fundamental 
theoretical  concepts  have  been  formulated  in  Hilbert  space  [6,7]. 

The  state  of  a  discretized  field  can  be  represented  by  a  vector  in  the  Hilbert 
space.  The  specification  of  the  mesh  node  connections  and  the  boundary  con- 
'Lehrstuhl  fur  Hochfcequenztechnik,  Teebnisebe  Univecsitat  Munchen,  Ar- 
cisstrasse  21,  D-8000  Munich  2,  Fed.  Rep.  Germany 
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ditions  is  done  by  operetois  in  the  Hilbert  space.  The  Hilbert  space  representa¬ 
tion  also  allows  the  description  of  geometric  structures  by  projection  operators. 
The  space  and  time  evolution  of  the  field  state  vector  is  governed  by  operator 
equations. 

In  field  theory,  the  field  propagation  in  spatial  domains  may  be  treated  using 
Green’s  functions  [8].  The  concept  of  Green’s  functions  may  also  be  applied 
to  discrete  time  dommn  field  computation  [9].  Discrete  time  domain  Green’s 
functions  allow  to  model  the  relation  between  the  field  values  on  the  boundaries 
if  the  knowledge  of  the  field  in  the  spatial  domains  beyond  the  boundaries  is 
not  required. 

In  this  paper,  the  algebraic  computation  of  the  discrete  Green’s  function  is 
investigated.  Our  approach  is  based  on  a  Hilbert  space  representation  of  the 
space-  and  time  discretired  electromagnetic  field.  The  discrete  Green’s  func¬ 
tions  may  be  represented  by  a  Neumann  series  in  space-  and  time-shift  oper¬ 
ators.  The  system  of  difference  equations  governing  the  time  evolution  of  the 
electromagnetic  field  in  configuration  space  in  derived  from  the  operator  equa¬ 
tion  for  the  field  state  vector  in  the  Hilbert  lipace.  First  results  are  presented 
for  the  two-dimrasional  case. 


2  The  Two-dimensional  TLM  Method 

The  electromagnetic  field  is  discretized  within  space  and  time.  The  space 
is  modelled  by  a  mesh  of  tran$mis»on  lines  connecting  the  sample  points  in 
space.  The  held  computation  algorithm  consists  of  two  steps; 

•  The  propagation  of  wave  pulses  fmm  the  mesh  nodes  to  the  neighbouring 
nodes. 

•  The  scattering  of  the  wave  puls^  in  the  mesh  nodes. 

In  the  following,  we  restrict  our  considerations  to  the  two^imensional  case 
with  the  transverse  electric  held.  In  the  shunt  TLM  model,  voltage  wave 
amplitudes  are  used  instead  of  total  voltage  and  current.  The  voltage  wave 
amplitudes  of  the  incident  and  the  reflected  waves  are  given  by  k<im.n  ^d 
khm,n.  The  left  index  k  denotes  the  discrete  time  coordinate  and  the  right 
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indices  m  &nd  n  denote  the  two  discrete  sp&ce  coordinates.  We  consider  the 
TLM  mesh  to  be  composed  by  elementary  TLM  shunt  node  fou]>port$  as 
shown  in  Fig.  1,  where  each  of  the  four  arms  is  of  length  A//2.  The  scattering 
in  this  elementary  four-port  is  connected  with  the  time  delay  At. 

The  scattering  of  the  wave  pulses  is  described  by 


with  the  scattering  matrix  5  given  by 

f-l  I  1  11 

1  .1  1  1 

5=  I  I  (2) 

2  3  2  3 

If  }  \  -U 

With  the  icattenng,  s  time  delay  of  At  is  usocisted  and  therefore,  the  time 
index  k  is  incremented  by  one.  The  scattered  pulses  are  the  incident  pulses  of 
the  neighbouring  elementary  cell.  This  is  described  by 

SUs.m.n  ~  *^^.,.,.+1 


3  The  Discrete  Field  State  Space 

In  the  TLM  model,  the  field  state  at  a  given  discrete  time  is  described  com¬ 
pletely  by  specifying  the  amplitudes  of  the  lout  wave  pulses  incident  to  each 
mesh  node .  The  space  of  the  voltage  wave  amplitudes  of  the  incident  and  the 
reflected  waves  and  is  the  fou^dimensional  real  vector  space 

R*.  In  order  to  develop  our  formalism  in  a  mote  general  way  we  introduce  the 
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Figure  li  A  two-dimensional  TLM  shunt  node  four-port. 

four-dimensional  complex  vector  space  C*  for  representing  the  wave  amplitudes 
SOn.n  sod 

In  order  to  describe  the  whole  mesh  state,  we  introduce  the  Hilbert  space  'Hm 
which  allows  to  map  each  mesh  node  onto  an  ortonormal  set  of  base  vectors  of 
The  time  states  are  represented  by  the  Hilbert  space  W|.  With  each  pair 
of  discrete  spatiai  coordinates  (m.ri),  a  basis  vector  of  tin  is  associated  and 
with  each  k,  a  basis  vector  of  Hi  is  associated.  We  now  introduce  the  state 
space  H  given  by  the  Cartesian  product  of  C*,  H»  and  Hi 

H  =  C‘®H«®H,  (4) 

The  space  H  is  a  Hilbert  space,  too.  The  complete  time  evolution  of  the  field 
state  within  the  whole  three-dimensional  space-time  may  now  be  represented 
by  a  single  vector  in  H.  Using  the  bra-ket  notation  introduced  by  Dirac  (6), 
the  orthonormal  basis  vectors  of  H  are  given  by  the  bra-vectors  |hi  m,  n).  The 
ket-vector  {it;m,n|  is  the  Hermitian  conjugate  of  |h;m,n).  The  orthogonality 
relations  are  given  by 

{^Ijrni,nilh3;m2,nr)  w  (5) 
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The  incident  and  reflected  voltage  waves  ate  represented  by 

llfcim,n) 


+»  +flO  +00 

W=  E  E  E 

Jts«oo  ma^oo  fls^oo 


and 


+00  +00  +00 

E  E  E 

Aa^oo  ma^oona^oo 


n) 


(6) 


(7) 


in  the  Hilbert  space  H.  We  define  the  shift  operators  JT,  Y  and  their  Hermitian 
conjugates  X*  and  V'  by 


JC|itjm,n)  =  H;m  +  l,n) 
jr'liimin)  =  Hjm-l,n) 
y|ibjm,n)  =  |ib;m,n  +  l) 
y*litjm,n)  =  |tjm,n  — 1) 


(8) 


The  operators  X  and  y  shift  the  field  state  by  one  interval  Af  in  the  positive 
m-  and  n-direction,  respectively.  Their  Hermitian  conjugates  X*  and  y'  shift 
the  field  state  in  the  opposite  direction. 


We  define  the  time  shift  operator  T.  The  time  shift  operator  increments  k  by 
f  i.e.  it  shifts  the  field  state  by  At  in  the  positive  time  direction.  If  the  time 
shift  operator  is  applied  to  a  vector  {Jl;m,n),  we  obtain 

T  |i;m,n)  =  It  + 1; m.n)  (9) 


We  introduce  the  connection  operator  F  jfiven  by 

0X00 
X'  0  0  0 

“  0  0  0  y 

0  0  y  0 


(10) 


With  the  connection  operator  F,  equation  (3)  yields  the  operator  equation 

|h)  =  r|<t)  (11) 
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Figure  2:  A  sputial  domain  witiun  the  TLM  mesh. 

The  inner  domain  projection  operator  projects  the  drcuit  space  H  on  the  inner 
ports  o{  the  domain  D,  Since  the  projection  operator  P/  and  the  connection 
operator  P  are  commuting,  i.e. 

iP/,rl  =  0  (22) 

we  obtain 

W,  =  rW,  (23) 

Applying  <iiakoptic$  to  TLM  structures  requires  the  computetion  of  the  wave 
pulses  scattered  at  the  domain  boundaries.  The  initial  conditions  or  boundary 
conditions  are  given  by  the  wave  pulses  inddent  on  the  boundary  ports.  We 
apply  the  projection  operators  PfPo  and  PbPd  in  order  to  separate  the  field 
states  |a)  and  |&)  into  the  inner  field  states  (a),  and  |&)j  and  the  boundary  states 
\a)g  and  \h)^.  PVom  eq.  (14)  we  obtain 

\h)B  =  TSBB  la)B  +  TSfl,  |a), 

16),  =  r5,B  |a)B  +  T5/,  la), 

with 

Sbb  =  Pb  S  Pb 
Sbi  -PbS  Pj 
Sib  -PiSPb 
Sn^PiSPi 


(25) 


Figure  3:  The  inner  ports  of  a  TLM  domain. 

Using  eqs.  (23)  and  (24),  we  eliminate  the  inner  domain  states  |o);  and  |&); 
and  obtain 

|4)b  =  [rSfis  +  TSst  (1  -  rrs,,)-'  rrs/s)  la)^  (26) 

This  is  the  relation  between  the  incident  and  scattered  boundary  state.  It 
describes  the  evolution  of  the  boundary  field  state  without  knowledge  of  the 
inner  field  state.  It  has  to  be  considered  that  the  operator  equation  (26) 
is  nonlocal  with  respect  to  both  space  and  time.  We  expand  the  operator 
(1  -  rT5//)'*  into  a  Neumann  series  (10,11)  and  obtain 

(i-rTS„)-'  =  f;T'(r5„)'  (27) 

r»o 

Inserting  this  into  eq.  (26)  yields  the  boundary  staU  evolution  equation 

\h)B^G\a)s  (28) 

with  the  boundary  field  evolution  operator  G  ^ven  by 

G  =  |r5B8  +  Sb,  (L3^«(r5„)')  rs/aj  (29) 

The  boundary  field  operator  G  gives  the  relation  between  the  boundary  state 
vector  (a}f  representing  the  wave  pulses  incident  on  the  boundary  and  the 


boundary  state  vector  i6)£  representing  the  wave  pulses  reflected  through  the 
boundary.  Eq.  (28)  is  the  general  formulation  of  the  boundary  element  problem 
in  the  Hilbert  space.  Since  the  Neumann  series  is  an  inflnite  geometrical  series 
in  space*  and  time*shift  operators,  the  boundary  field  operator  is  nonlocal 
with  respect  to  space  and  time. 


4  The  Discrete  Two-dimensional  Green’s 
function 

As  an  example,  we  derive  the  discrete  Green's  function  for  the  half*plane.  The 
discrete  Green's  function  for  the  balf*plane  is  given  by  the  projection  of  the 
boundary  state  evolution  operator  equation  (28)  onto  configuration  space  for 
a  point-like  initial  state  \a)g.  The  balf-pluie  (Fig.  4)  is  defined  by  the  domiun 
projection  operator  Pd  given  by 

Pd  =  Z  (30) 

k,n  maO 


As  in  the  shunt  TLM-model,  voltage  wave  amplitudes  instead  of  total  voltages 
are  used,  a  new  Green’s  function  for  wave  amplitudes  has  to  be  defined.  For 
a  boundary  problem,  the  Green's  function  in  discrete  formulation  is  given  by 
the  convolution 

k^u  =  kGn  *  *Oii'  (31) 


where  ik  Oa  is  the  column  vector  of  the  incident  impulse  functions  at  the  time 
k^t  and  at  the  boundary  node  number  n.  «  bn  is  the  column  vector  of  the 
scattered  output  wave  impulses  at  the  time  kAt  and  at  the  boundary 
node,  k  Gn  is  the  discrete  Green's  function  for  an  arbitrary  boundary  with 
n  boundary  nodes.  It  describes  the  relation  between  the  incident  anv  the 
scattered  wave  amplitudes  in  the  boundary  ports. 

For  the  half-plane,  the  boundary  is  given  by  m  =  0  and  n  = 
-00.:., ,— 1,0,1, ..,.00.  Thereforeeq,  (31)yields 


t  t  .-.'Gw 

fc'v-oo 


t 

t 

I 


e  <*n’ 


(32) 
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Figure  i:  The  homogeneous  two-dimensional  half-space. 

The  boundary  state  evolution  enuation  (28)  may  be  expressed  by  the  discrete 
Green’s  function,  eq,  (32),  via 


Wa  =  G|->8 


(33) 


where  the  boundary  held  evolution  operator  is  given  by 


10  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 


S  s-fG„-«>|i;0,n)(i'i0,n'| 


(34) 


In  order  to  calculate  the  Green's  function  for  the  boundary  of  the  half-plane, 
we  start  from  an  impulsive  excitation  at  n'  =  0,  IF  =  0  given  by 


(35) 


Mapping  eqs.  (28)  with  (29)  ard  (35),  (33)  with  (34)  and  (35)  to  configuration 
space  by  multiplying  both  equations  from  the  left  side  with  (i;  0,  n|,  we  obtrun 
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a  system  of  partial  difference  equations  which  can  be  solved  by  transforming 
it  to  frequency-  and  momentum-space. 


We  obt^un  an  algebraic  expression  for  the  Greenes  function 

kGn  =  2^*,n+l  —  2^a,n-I  +  ^  a+l/n  -  ^  jfc-s/ft 

*-i  1 


,ma  ®  * 


*-2  j  j 

j«0  ®  ^ 


+  I 

with  n  K  0, l,2,.x.oo;  Je  =  2, 3, 4,. ...oo  and 
*G-„=5*G« 

for  n  <  0. 

The  function  is  given  by 

*  (VJ)*.! 


(35) 


(37) 


irj\  2r  j((_2«  +  r-nj 

Id  Fij.  4,  ^  G„  is  depicted  for  n  =  -9..., -1,0,1,... 9;  i  =  1,2,.;. .  10. 


(38) 
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Figure  5:  Values  of  the  Greenes  function 
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Solving  eigenvalue  and  steady-state  problems 
using  time-domain  models 


Gunnar  Nitsche 

Lelirstuhl  fur  Nachrichtentechnik 
Ruhr-Universitat  Bochum 
D-4630  Bochum  1 
Germeuiy 


Time-domain  modelling  of  partial  differential  equations  has  become 
populM  in  the  last  years  for  several  reasons.  The  detailed  time-domain 
behaviour,  however,  is  often  not  of  primary  concern,  but  one  is  more  inter¬ 
ested  in  the  eigenvalues  and  eigenfunctions  of  a  system  or  its  steady-state 
response  to  a  sinusoidal  excitation.  These  kinds  of  problems  are  usually 
approached  by  inethods  based  on  the  discrete  Fourier  transform. 

A  new  alternative  approach  to  efficiently  compute  the  low  frequency 
eigeumodes  of  a  time-domain  model  approximating  a  physical  system  will 
be  proposed.  The  method  is  based  on  principles  known  from  digital  signal 
processing,  in  particular  from  parametric  spectrum  estimation,  so  it  is  not 
surprising  that  the  achievable  accuracy  is  much  higher  than  the  accuracy 
of  the  uou-parametrie  Fourier  transform  approach.  The  algoritlim  works  in 
the^  general  lossy  case  even  for  a  very  large  number  of  unknowns  and  can 
easily  be  extended  to  calculate  steady-state  solutions  for  several  different 
frequencies  simultaneously. 


Restaurant  Guide 

Near  the  TU  there  are  some  restaurants  and  coffee-houses  were  you  can  go  to.  Some 
of  these  are  listed  below.  The  number  is  rdated  to  the  numbers  on  the  map.  So  you 
can  easily  find  your  location. 

1.  CANTON  chtntst  restaurant,  Theresienstr.  49 

2.  BEI  MARIO  pizzeria,  italian  restaurant,  Lmsenstr. 

3.  BELLA  ITALIA  pizzerta,  italian  rtsiaurant,  Turkenstr. 

4.  HIEN  Vietnamese  food,  Schellingstr.  91 

5.  CAFE  ALTSCHWABING  coffee  house,  bistro,  Schellingstr. 

6.  WEINSCHATULLE  restaurant,  Thei^ienstr. 

7.  MC  DONALDS  fast  food,  Augustenstr. 


On  thursday  evening  there  is  a  sodai  event  consisting  of  a  Song  Recital  with  Piano 
and  dinner  for  all. workshop  participants  at  the  restaurant  “Seehaus”,  Englischer 
Garten  (see  the  map  below).  We  recommend  you  to  reach  the  Seehaus  by  car  or  taxi 
coming  &om  the  Mittlerer  Ring. 


